Changes in / [85da4e:f42e95]


Ignore:
Location:
src
Files:
4 added
25 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    r85da4e rf42e95  
    88ANALYSISHEADER = analysis_bonds.hpp analysis_correlation.hpp
    99
    10 SOURCE = ${ANALYSISSOURCE} ${ATOMSOURCE} bond.cpp bondgraph.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp errorlogger.cpp graph.cpp helpers.cpp info.cpp leastsquaremin.cpp linkedcell.cpp log.cpp logger.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp vector.cpp verbose.cpp
    11 HEADER = ${ANALYSISHEADER} ${ATOMHEADER} bond.hpp bondgraph.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp errorlogger.hpp graph.hpp helpers.hpp info.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp log.hpp logger.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp vector.hpp verbose.hpp
     10SOURCE = ${ANALYSISSOURCE} ${ATOMSOURCE} bond.cpp bondgraph.cpp boundary.cpp config.cpp element.cpp ellipsoid.cpp errorlogger.cpp graph.cpp helpers.cpp info.cpp leastsquaremin.cpp linkedcell.cpp log.cpp logger.cpp memoryusageobserver.cpp moleculelist.cpp molecule.cpp molecule_dynamics.cpp molecule_fragmentation.cpp molecule_geometry.cpp molecule_graph.cpp molecule_pointcloud.cpp parser.cpp periodentafel.cpp tesselation.cpp tesselationhelpers.cpp triangleintersectionlist.cpp vector.cpp verbose.cpp World.cpp
     11HEADER = ${ANALYSISHEADER} ${ATOMHEADER} bond.hpp bondgraph.hpp boundary.hpp config.hpp defs.hpp element.hpp ellipsoid.hpp errorlogger.hpp graph.hpp helpers.hpp info.hpp leastsquaremin.hpp linkedcell.hpp lists.hpp log.hpp logger.hpp memoryallocator.hpp memoryusageobserver.hpp molecule.hpp molecule_template.hpp parser.hpp periodentafel.hpp stackclass.hpp tesselation.hpp tesselationhelpers.hpp triangleintersectionlist.cpp vector.hpp verbose.hpp World.hpp
    1212
    1313BOOST_LIB = $(BOOST_LDFLAGS) $(BOOST_MPL_LIB)

  • src/analysis_correlation.cpp

    r85da4e rf42e95  
    1515#include "tesselation.hpp"
    1616#include "tesselationhelpers.hpp"
     17#include "triangleintersectionlist.hpp"
    1718#include "vector.hpp"
    1819#include "verbose.hpp"
     20#include "World.hpp"
    1921
    2022
    … …  
    5557                if (Walker->nr < OtherWalker->nr)
    5658                  if ((type2 == NULL) || (OtherWalker->type == type2)) {
    57                     distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
     59                    distance = Walker->node->PeriodicDistance(OtherWalker->node, World::get()->cell_size);
    5860                    //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
    5961                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
    … …  
    9698  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    9799    if ((*MolWalker)->ActiveFlag) {
    98       double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
     100      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
    99101      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    100102      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    … …  
    174176        Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
    175177        if ((type == NULL) || (Walker->type == type)) {
    176           distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
     178          distance = Walker->node->PeriodicDistance(point, World::get()->cell_size);
    177179          Log() << Verbose(4) << "Current distance is " << distance << "." << endl;
    178180          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
    … …  
    208210  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    209211    if ((*MolWalker)->ActiveFlag) {
    210       double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
     212      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
    211213      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    212214      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    … …  
    255257
    256258  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
    257     Log() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
     259    eLog() << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
    258260    return outmap;
    259261  }
    … …  
    261263  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    262264    if ((*MolWalker)->ActiveFlag) {
    263       Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    264       atom *Walker = (*MolWalker)->start;
    265       while (Walker->next != (*MolWalker)->end) {
    266         Walker = Walker->next;
    267         Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl;
     265      Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl;
     266      atom *Walker = (*MolWalker)->start;
     267      while (Walker->next != (*MolWalker)->end) {
     268        Walker = Walker->next;
     269        //Log() << Verbose(1) << "Current atom is " << *Walker << "." << endl;
    268270        if ((type == NULL) || (Walker->type == type)) {
    269           triangle = Surface->FindClosestTriangleToVector(Walker->node, LC );
    270           if (triangle != NULL) {
    271             distance = DistanceToTrianglePlane(Walker->node, triangle);
    272             outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
    273           }
    274         }
    275       }
    276     }
     271          TriangleIntersectionList Intersections(Walker->node,Surface,LC);
     272          distance = Intersections.GetSmallestDistance();
     273          triangle = Intersections.GetClosestTriangle();
     274          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
     275        }
     276      }
     277    } else
     278      Log() << Verbose(1) << "molecule " << (*MolWalker)->name << " is not active." << endl;
     279
    277280
    278281  return outmap;
    … …  
    312315  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    313316    if ((*MolWalker)->ActiveFlag) {
    314       double * FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
     317      double * FullMatrix = ReturnFullMatrixforSymmetric(World::get()->cell_size);
    315318      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    316319      Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    … …  
    372375  *file << "BinStart\tCount" << endl;
    373376  for (BinPairMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    374     *file << runner->first << "\t" << runner->second << endl;
     377    *file << setprecision(8) << runner->first << "\t" << runner->second << endl;
    375378  }
    376379};
    … …  
    385388  *file << "BinStart\tAtom1\tAtom2" << endl;
    386389  for (PairCorrelationMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    387     *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
     390    *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
    388391  }
    389392};
    … …  
    400403    *file << runner->first;
    401404    for (int i=0;i<NDIM;i++)
    402       *file << "\t" << (runner->second.first->node->x[i] - runner->second.second->x[i]);
     405      *file << "\t" << setprecision(8) << (runner->second.first->node->x[i] - runner->second.second->x[i]);
    403406    *file << endl;
    404407  }
    … …  
    413416  Info FunctionInfo(__func__);
    414417  *file << "BinStart\tTriangle" << endl;
    415   for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
    416     *file << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
    417   }
    418 };
    419 
     418  if (!map->empty())
     419    for (CorrelationToSurfaceMap::const_iterator runner = map->begin(); runner != map->end(); ++runner) {
     420      *file << setprecision(8) << runner->first << "\t" << *(runner->second.first) << "\t" << *(runner->second.second) << endl;
     421    }
     422};
     423

  • src/atom_bondedparticle.cpp

    r85da4e rf42e95  
    5656 * \param *AdjacencyFile output stream
    5757 */
    58 void BondedParticle::OutputAdjacency(ofstream *AdjacencyFile) const
     58void BondedParticle::OutputAdjacency(ofstream * const AdjacencyFile) const
    5959{
    6060  *AdjacencyFile << nr << "\t";
    … …  
    6262    *AdjacencyFile << (*Runner)->GetOtherAtom(this)->nr << "\t";
    6363  *AdjacencyFile << endl;
     64};
     65
     66/** Output of atom::nr along each bond partner per line.
     67 * Only bonds are printed where atom::nr is smaller than the one of the bond partner.
     68 * \param *AdjacencyFile output stream
     69 */
     70void BondedParticle::OutputBonds(ofstream * const BondFile) const
     71{
     72  for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner))
     73    if (nr < (*Runner)->GetOtherAtom(this)->nr)
     74      *BondFile << nr << "\t" << (*Runner)->GetOtherAtom(this)->nr << "\n";
    6475};
    6576

  • src/atom_bondedparticle.hpp

    r85da4e rf42e95  
    4444  int CorrectBondDegree();
    4545  void OutputBondOfAtom() const;
    46   void OutputAdjacency(ofstream *AdjacencyFile) const;
     46  void OutputAdjacency(ofstream * const AdjacencyFile) const;
     47  void OutputBonds(ofstream * const BondFile) const;
    4748  void OutputOrder(ofstream *file) const;
    4849

  • src/boundary.cpp

    r85da4e rf42e95  
    1717#include "tesselation.hpp"
    1818#include "tesselationhelpers.hpp"
     19#include "World.hpp"
    1920
    2021#include<gsl/gsl_poly.h>
    … …  
    789790 * \param *filler molecule which the box is to be filled with
    790791 * \param configuration contains box dimensions
     792 * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
    791793 * \param distance[NDIM] distance between filling molecules in each direction
    792794 * \param boundary length of boundary zone between molecule and filling mollecules
    … …  
    797799 * \return *mol pointer to new molecule with filled atoms
    798800 */
    799 molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
     801molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
    800802{
    801803        Info FunctionInfo(__func__);
    … …  
    804806  int N[NDIM];
    805807  int n[NDIM];
    806   double *M =  ReturnFullMatrixforSymmetric(filler->cell_size);
     808  double *M =  ReturnFullMatrixforSymmetric(World::get()->cell_size);
    807809  double Rotations[NDIM*NDIM];
     810  double *MInverse = InverseMatrix(M);
    808811  Vector AtomTranslations;
    809812  Vector FillerTranslations;
    810813  Vector FillerDistance;
     814  Vector Inserter;
    811815  double FillIt = false;
    812816  atom *Walker = NULL;
    813817  bond *Binder = NULL;
    814   int i = 0;
    815   LinkedCell *LCList[List->ListOfMolecules.size()];
    816818  double phi[NDIM];
    817   class Tesselation *TesselStruct[List->ListOfMolecules.size()];
    818 
    819   i=0;
    820   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    821     Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
    822     LCList[i] = new LinkedCell((*ListRunner), 10.); // get linked cell list
    823     Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
    824     TesselStruct[i] = NULL;
    825     FindNonConvexBorder((*ListRunner), TesselStruct[i], (const LinkedCell *&)LCList[i], 5., NULL);
    826     i++;
    827   }
     819  map<molecule *, Tesselation *> TesselStruct;
     820  map<molecule *, LinkedCell *> LCList;
     821
     822  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
     823    if ((*ListRunner)->AtomCount > 0) {
     824      Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl;
     825      LCList[(*ListRunner)] = new LinkedCell((*ListRunner), 10.); // get linked cell list
     826      Log() << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
     827      TesselStruct[(*ListRunner)] = NULL;
     828      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell *&)LCList[(*ListRunner)], 5., NULL);
     829    }
    828830
    829831  // Center filler at origin
    830   filler->CenterOrigin();
     832  filler->CenterEdge(&Inserter);
    831833  filler->Center.Zero();
     834  Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl;
     835  Binder = filler->first;
     836  while(Binder->next != filler->last) {
     837    Binder = Binder->next;
     838    Log() << Verbose(2) << "  " << *Binder << endl;
     839  }
    832840
    833841  filler->CountAtoms();
    … …  
    851859        CurrentPosition.Init((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    852860        CurrentPosition.MatrixMultiplication(M);
    853         Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "." << endl;
    854         // Check whether point is in- or outside
    855         FillIt = true;
    856         i=0;
    857         for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    858           // get linked cell list
    859           if (TesselStruct[i] == NULL) {
    860             eLog() << Verbose(0) << "TesselStruct of " << (*ListRunner) << " is NULL. Didn't we pre-create it?" << endl;
    861             FillIt = false;
     861        // create molecule random translation vector ...
     862        for (int i=0;i<NDIM;i++)
     863          FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
     864        Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << "." << endl;
     865
     866        // go through all atoms
     867        for (int i=0;i<filler->AtomCount;i++)
     868          CopyAtoms[i] = NULL;
     869        Walker = filler->start;
     870        while (Walker->next != filler->end) {
     871          Walker = Walker->next;
     872
     873          // create atomic random translation vector ...
     874          for (int i=0;i<NDIM;i++)
     875            AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
     876
     877          // ... and rotation matrix
     878          if (DoRandomRotation) {
     879            for (int i=0;i<NDIM;i++) {
     880              phi[i] = rand()/(RAND_MAX/(2.*M_PI));
     881            }
     882
     883            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
     884            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
     885            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
     886            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
     887            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
     888            Rotations[7] =               sin(phi[1])                                                ;
     889            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
     890            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
     891            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
     892          }
     893
     894          // ... and put at new position
     895          Inserter.CopyVector(&(Walker->x));
     896          if (DoRandomRotation)
     897            Inserter.MatrixMultiplication(Rotations);
     898          Inserter.AddVector(&AtomTranslations);
     899          Inserter.AddVector(&FillerTranslations);
     900          Inserter.AddVector(&CurrentPosition);
     901
     902          // check whether inserter is inside box
     903          Inserter.MatrixMultiplication(MInverse);
     904          FillIt = true;
     905          for (int i=0;i<NDIM;i++)
     906            FillIt = FillIt && (Inserter.x[i] >= -MYEPSILON) && ((Inserter.x[i]-1.) <= MYEPSILON);
     907          Inserter.MatrixMultiplication(M);
     908
     909          // Check whether point is in- or outside
     910          for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
     911            // get linked cell list
     912            if (TesselStruct[(*ListRunner)] != NULL) {
     913              const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
     914              FillIt = FillIt && (distance > boundary) && ((MaxDistance < 0) || (MaxDistance > distance));
     915            }
     916          }
     917          // insert into Filling
     918          if (FillIt) {
     919            Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl;
     920            // copy atom ...
     921            CopyAtoms[Walker->nr] = new atom(Walker);
     922            CopyAtoms[Walker->nr]->x.CopyVector(&Inserter);
     923            Filling->AddAtom(CopyAtoms[Walker->nr]);
     924            Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
    862925          } else {
    863             const double distance = (TesselStruct[i]->GetDistanceSquaredToSurface(CurrentPosition, LCList[i]));
    864             FillIt = FillIt && (distance > boundary*boundary);
    865             if (FillIt) {
    866               Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is outer point." << endl;
    867             } else {
    868               Log() << Verbose(1) << "INFO: Position at " << CurrentPosition << " is inner point or within boundary." << endl;
    869               break;
    870             }
    871             i++;
     926            Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl;
     927            CopyAtoms[Walker->nr] = NULL;
     928            continue;
    872929          }
    873930        }
    874 
    875         if (FillIt) {
    876           // fill in Filler
    877           Log() << Verbose(2) << "Space at " << CurrentPosition << " is unoccupied by any molecule, filling in." << endl;
    878 
    879           // create molecule random translation vector ...
    880           for (int i=0;i<NDIM;i++)
    881             FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
    882           Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;
    883 
    884           // go through all atoms
    885           Walker = filler->start;
    886           while (Walker->next != filler->end) {
    887             Walker = Walker->next;
    888             // copy atom ...
    889             CopyAtoms[Walker->nr] = new atom(Walker);
    890 
    891             // create atomic random translation vector ...
    892             for (int i=0;i<NDIM;i++)
    893               AtomTranslations.x[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
    894 
    895             // ... and rotation matrix
    896             if (DoRandomRotation) {
    897               for (int i=0;i<NDIM;i++) {
    898                 phi[i] = rand()/(RAND_MAX/(2.*M_PI));
    899               }
    900 
    901               Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
    902               Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
    903               Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
    904               Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
    905               Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
    906               Rotations[7] =               sin(phi[1])                                                ;
    907               Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
    908               Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
    909               Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
    910             }
    911 
    912             // ... and put at new position
    913             if (DoRandomRotation)
    914               CopyAtoms[Walker->nr]->x.MatrixMultiplication(Rotations);
    915             CopyAtoms[Walker->nr]->x.AddVector(&AtomTranslations);
    916             CopyAtoms[Walker->nr]->x.AddVector(&FillerTranslations);
    917             CopyAtoms[Walker->nr]->x.AddVector(&CurrentPosition);
    918 
    919             // insert into Filling
    920 
    921             // FIXME: gives completely different results if CopyAtoms[..] used instead of Walker, why???
    922             Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl;
    923             Filling->AddAtom(CopyAtoms[Walker->nr]);
    924           }
    925 
    926           // go through all bonds and add as well
    927           Binder = filler->first;
    928           while(Binder->next != filler->last) {
    929             Binder = Binder->next;
    930             Log() << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
     931        // go through all bonds and add as well
     932        Binder = filler->first;
     933        while(Binder->next != filler->last) {
     934          Binder = Binder->next;
     935          if ((CopyAtoms[Binder->leftatom->nr] != NULL) && (CopyAtoms[Binder->rightatom->nr] != NULL)) {
     936            Log()  << Verbose(3) << "Adding Bond between " << *CopyAtoms[Binder->leftatom->nr] << " and " << *CopyAtoms[Binder->rightatom->nr]<< "." << endl;
    931937            Filling->AddBond(CopyAtoms[Binder->leftatom->nr], CopyAtoms[Binder->rightatom->nr], Binder->BondDegree);
    932938          }
    933         } else {
    934           // leave empty
    935           Log() << Verbose(2) << "Space at " << CurrentPosition << " is occupied." << endl;
    936939        }
    937940      }
    938941  Free(&M);
    939 
    940   // output to file
    941   TesselStruct[0]->LastTriangle = NULL;
    942   StoreTrianglesinFile(Filling, TesselStruct[0], "Tesselated", ".dat");
    943 
    944   for (size_t i=0;i<List->ListOfMolecules.size();i++) {
    945         delete(LCList[i]);
    946         delete(TesselStruct[i]);
    947   }
     942  Free(&MInverse);
     943
    948944  return Filling;
    949945};

  • src/boundary.hpp

    r85da4e rf42e95  
    4949
    5050double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
    51 molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
     51molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
    5252void FindConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
    5353Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);

  • src/builder.cpp

    r85da4e rf42e95  
    6868#include "periodentafel.hpp"
    6969#include "version.h"
     70#include "World.hpp"
    7071
    7172/********************************************* Subsubmenu routine ************************************/
    … …  
    103104        Log() << Verbose(0) << "Enter absolute coordinates." << endl;
    104105        first = new atom;
    105         first->x.AskPosition(mol->cell_size, false);
     106        first->x.AskPosition(World::get()->cell_size, false);
    106107        first->type = periode->AskElement();  // give type
    107108        mol->AddAtom(first);  // add to molecule
    … …  
    114115          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
    115116          Log() << Verbose(0) << "Enter reference coordinates." << endl;
    116           x.AskPosition(mol->cell_size, true);
     117          x.AskPosition(World::get()->cell_size, true);
    117118          Log() << Verbose(0) << "Enter relative coordinates." << endl;
    118           first->x.AskPosition(mol->cell_size, false);
     119          first->x.AskPosition(World::get()->cell_size, false);
    119120          first->x.AddVector((const Vector *)&x);
    120121          Log() << Verbose(0) << "\n";
    … …  
    131132          second = mol->AskAtom("Enter atom number: ");
    132133          Log() << Verbose(0) << "Enter relative coordinates." << endl;
    133           first->x.AskPosition(mol->cell_size, false);
     134          first->x.AskPosition(World::get()->cell_size, false);
    134135          for (int i=NDIM;i--;) {
    135136            first->x.x[i] += second->x.x[i];
    … …  
    358359    case 'b': // normal vector of mirror plane
    359360      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
    360       n.AskPosition(mol->cell_size,0);
     361      n.AskPosition(World::get()->cell_size,0);
    361362      n.Normalize();
    362363      break;
    … …  
    425426    case 'b': // normal vector of mirror plane
    426427      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
    427       n.AskPosition(mol->cell_size,0);
     428      n.AskPosition(World::get()->cell_size,0);
    428429      n.Normalize();
    429430      break;
    … …  
    655656static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
    656657{
    657   atom *first, *second;
     658  atom *first, *second, *third;
    658659  molecule *mol = NULL;
    659660  Vector x,y,z,n; // coordinates for absolute point in cell volume
    … …  
    667668  Log() << Verbose(0) << "r - remove an atom" << endl;
    668669  Log() << Verbose(0) << "b - scale a bond between atoms" << endl;
     670  Log() << Verbose(0) << "t - turn an atom round another bond" << endl;
    669671  Log() << Verbose(0) << "u - change an atoms element" << endl;
    670672  Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
    … …  
    748750      break;
    749751
     752    case 't': // turn/rotate atom
     753      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     754        if ((*ListRunner)->ActiveFlag) {
     755          mol = *ListRunner;
     756          Log() << Verbose(0) << "Turning atom around another bond - first is atom to turn, second (central) and third specify bond" << endl;
     757          first = mol->AskAtom("Enter turning atom: ");
     758          second = mol->AskAtom("Enter central atom: ");
     759          third  = mol->AskAtom("Enter bond atom: ");
     760          cout << Verbose(0) << "Enter new angle in degrees: ";
     761          double tmp = 0.;
     762          cin >> tmp;
     763          // calculate old angle
     764          x.CopyVector((const Vector *)&first->x);
     765          x.SubtractVector((const Vector *)&second->x);
     766          y.CopyVector((const Vector *)&third->x);
     767          y.SubtractVector((const Vector *)&second->x);
     768          double alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
     769          cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
     770          cout << Verbose(0) << alpha << " degrees" << endl;
     771          // rotate
     772          z.MakeNormalVector(&x,&y);
     773          x.RotateVector(&z,(alpha-tmp)*M_PI/180.);
     774          x.AddVector(&second->x);
     775          first->x.CopyVector(&x);
     776          // check new angle
     777          x.CopyVector((const Vector *)&first->x);
     778          x.SubtractVector((const Vector *)&second->x);
     779          alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
     780          cout << Verbose(0) << "new Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
     781          cout << Verbose(0) << alpha << " degrees" << endl;
     782        }
     783      break;
     784
    750785    case 'u': // change an atom's element
    751786      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    … …  
    831866          x.Zero();
    832867          y.Zero();
    833           y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
     868          y.x[abs(axis)-1] = World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
    834869          for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
    835870            x.AddVector(&y); // per factor one cell width further
    … …  
    854889            mol->Translate(&x);
    855890          }
    856           mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     891          World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    857892        }
    858893      }
    … …  
    911946        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
    912947        Log() << Verbose(0) << "Enter translation vector." << endl;
    913         x.AskPosition(mol->cell_size,0);
     948        x.AskPosition(World::get()->cell_size,0);
    914949        mol->Center.AddVector((const Vector *)&x);
    915950     }
    … …  
    9711006        // center at set box dimensions
    9721007        mol->CenterEdge(&center);
    973         mol->cell_size[0] = center.x[0];
    974         mol->cell_size[1] = 0;
    975         mol->cell_size[2] = center.x[1];
    976         mol->cell_size[3] = 0;
    977         mol->cell_size[4] = 0;
    978         mol->cell_size[5] = center.x[2];
     1008        double * const cell_size = World::get()->cell_size;
     1009        cell_size[0] = center.x[0];
     1010        cell_size[1] = 0;
     1011        cell_size[2] = center.x[1];
     1012        cell_size[3] = 0;
     1013        cell_size[4] = 0;
     1014        cell_size[5] = center.x[2];
    9791015        molecules->insert(mol);
    9801016      }
    … …  
    13871423  enum ConfigStatus configPresent = absent;
    13881424  clock_t start,end;
     1425  double MaxDistance = -1;
    13891426  int argptr;
    13901427  molecule *mol = NULL;
    … …  
    14121449            Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
    14131450            Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
    1414             Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl;
     1451            Log() << Verbose(0) << "\t-C <type> [params] <output> <bin output> <BinWidth> <BinStart> <BinEnd>\tPair Correlation analysis." << endl;
    14151452            Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    14161453            Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
    … …  
    14181455            Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
    14191456            Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
    1420             Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
     1457            Log() << Verbose(0) << "\t-F <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
     1458            Log() << Verbose(0) << "\t-FF <MaxDistance> <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
    14211459            Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
    14221460            Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
    14231461            Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl;
     1462            Log() << Verbose(0) << "\t-j\t<path> Store all bonds to file." << endl;
     1463            Log() << Verbose(0) << "\t-J\t<path> Store adjacency per atom to file." << endl;
    14241464            Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
    14251465            Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    … …  
    14401480            Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl;
    14411481            Log() << Verbose(0) << "\t-V\t\tGives version information." << endl;
     1482            Log() << Verbose(0) << "\t-X\t\tset default name of a molecule." << endl;
    14421483            Log() << Verbose(0) << "Note: config files must not begin with '-' !" << endl;
    14431484            return (1);
    … …  
    14781519            Log() << Verbose(0) << "I won't parse trajectories." << endl;
    14791520            configuration.FastParsing = true;
     1521            break;
     1522          case 'X':
     1523            {
     1524              char **name = &(World::get()->DefaultName);
     1525              delete[](*name);
     1526              const int length = strlen(argv[argptr]);
     1527              *name = new char[length+2];
     1528              strncpy(*name, argv[argptr], length);
     1529              Log() << Verbose(0) << "Default name of new molecules set to " << *name << "." << endl;
     1530            }
    14801531            break;
    14811532          default:   // no match? Step on
    … …  
    16691720                  }
    16701721              }
    1671               if (mol == NULL) {
     1722              if ((mol == NULL) && (!molecules->ListOfMolecules.empty())) {
    16721723                mol = *(molecules->ListOfMolecules.begin());
    1673                 mol->ActiveFlag = true;
     1724                if (mol != NULL)
     1725                  mol->ActiveFlag = true;
    16741726              }
    16751727              break;
    16761728            case 'C':
    16771729              if (ExitFlag == 0) ExitFlag = 1;
    1678               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-') || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-')) {
     1730              if ((argptr >= argc)) {
    16791731                ExitFlag = 255;
    1680                 eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <Z> <output> <bin output>" << endl;
     1732                eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output> <BinStart> <BinEnd>" << endl;
    16811733                performCriticalExit();
    16821734              } else {
    1683                 ofstream output(argv[argptr+1]);
    1684                 ofstream binoutput(argv[argptr+2]);
    1685                 const double radius = 5.;
    1686 
    1687                 // get the boundary
    1688                 class molecule *Boundary = NULL;
    1689                 class Tesselation *TesselStruct = NULL;
    1690                 const LinkedCell *LCList = NULL;
    1691                 // find biggest molecule
    1692                 int counter  = 0;
    1693                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    1694                   if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
    1695                     Boundary = *BigFinder;
    1696                   }
    1697                   counter++;
     1735                switch(argv[argptr][0]) {
     1736                  case 'E':
     1737                    {
     1738                      if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+2])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-')) {
     1739                        ExitFlag = 255;
     1740                        eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C E <Z1> <Z2> <output> <bin output>" << endl;
     1741                        performCriticalExit();
     1742                      } else {
     1743                        ofstream output(argv[argptr+3]);
     1744                        ofstream binoutput(argv[argptr+4]);
     1745                        const double BinStart = atof(argv[argptr+5]);
     1746                        const double BinEnd = atof(argv[argptr+6]);
     1747
     1748                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
     1749                        element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+2]));
     1750                        PairCorrelationMap *correlationmap = PairCorrelation(molecules, elemental, elemental2);
     1751                        //OutputCorrelationToSurface(&output, correlationmap);
     1752                        BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
     1753                        OutputCorrelation ( &binoutput, binmap );
     1754                        output.close();
     1755                        binoutput.close();
     1756                        delete(binmap);
     1757                        delete(correlationmap);
     1758                        argptr+=7;
     1759                      }
     1760                    }
     1761                    break;
     1762
     1763                  case 'P':
     1764                    {
     1765                      if ((argptr+8 >= argc) || (!IsValidNumber(argv[argptr+1])) ||  (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+7])) || (!IsValidNumber(argv[argptr+8])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-') || (argv[argptr+5][0] == '-') || (argv[argptr+6][0] == '-')) {
     1766                        ExitFlag = 255;
     1767                        eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C P <Z1> <x> <y> <z> <output> <bin output>" << endl;
     1768                        performCriticalExit();
     1769                      } else {
     1770                        ofstream output(argv[argptr+5]);
     1771                        ofstream binoutput(argv[argptr+6]);
     1772                        const double BinStart = atof(argv[argptr+7]);
     1773                        const double BinEnd = atof(argv[argptr+8]);
     1774
     1775                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
     1776                        Vector *Point = new Vector((const double) atof(argv[argptr+1]),(const double) atof(argv[argptr+2]),(const double) atof(argv[argptr+3]));
     1777                        CorrelationToPointMap *correlationmap = CorrelationToPoint(molecules, elemental, Point);
     1778                        //OutputCorrelationToSurface(&output, correlationmap);
     1779                        BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
     1780                        OutputCorrelation ( &binoutput, binmap );
     1781                        output.close();
     1782                        binoutput.close();
     1783                        delete(Point);
     1784                        delete(binmap);
     1785                        delete(correlationmap);
     1786                        argptr+=9;
     1787                      }
     1788                    }
     1789                    break;
     1790
     1791                  case 'S':
     1792                    {
     1793                      if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-')) {
     1794                        ExitFlag = 255;
     1795                        eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C S <Z> <output> <bin output> <BinWidth> <BinStart> <BinEnd>" << endl;
     1796                        performCriticalExit();
     1797                      } else {
     1798                        ofstream output(argv[argptr+2]);
     1799                        ofstream binoutput(argv[argptr+3]);
     1800                        const double radius = 4.;
     1801                        const double BinWidth = atof(argv[argptr+4]);
     1802                        const double BinStart = atof(argv[argptr+5]);
     1803                        const double BinEnd = atof(argv[argptr+6]);
     1804                        double LCWidth = 20.;
     1805                        if (BinEnd > 0) {
     1806                          if (BinEnd > 2.*radius)
     1807                              LCWidth = BinEnd;
     1808                          else
     1809                            LCWidth = 2.*radius;
     1810                        }
     1811
     1812                        // get the boundary
     1813                        class molecule *Boundary = NULL;
     1814                        class Tesselation *TesselStruct = NULL;
     1815                        const LinkedCell *LCList = NULL;
     1816                        // find biggest molecule
     1817                        int counter  = 0;
     1818                        for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
     1819                          if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
     1820                            Boundary = *BigFinder;
     1821                          }
     1822                          counter++;
     1823                        }
     1824                        bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
     1825                        counter = 0;
     1826                        for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
     1827                          Actives[counter++] = (*BigFinder)->ActiveFlag;
     1828                          (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
     1829                        }
     1830                        LCList = new LinkedCell(Boundary, LCWidth);
     1831                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
     1832                        FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
     1833                        //int ranges[NDIM] = {1,1,1};
     1834                        CorrelationToSurfaceMap *surfacemap = CorrelationToSurface( molecules, elemental, TesselStruct, LCList); // for Periodic..(): ..., ranges );
     1835                        OutputCorrelationToSurface(&output, surfacemap);
     1836                        // check whether radius was appropriate
     1837                        {
     1838                        double start; double end;
     1839                        GetMinMax( surfacemap, start, end);
     1840                        if (LCWidth < end)
     1841                          eLog() << Verbose(1) << "Linked Cell width is smaller than the found range of values! Bins can only be correct up to: " << radius << "." << endl;
     1842                        }
     1843                        BinPairMap *binmap = BinData( surfacemap, BinWidth, BinStart, BinEnd );
     1844                        OutputCorrelation ( &binoutput, binmap );
     1845                        output.close();
     1846                        binoutput.close();
     1847                        for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
     1848                          (*BigFinder)->ActiveFlag = Actives[counter++];
     1849                        Free(&Actives);
     1850                        delete(LCList);
     1851                        delete(TesselStruct);
     1852                        delete(binmap);
     1853                        delete(surfacemap);
     1854                        argptr+=7;
     1855                      }
     1856                    }
     1857                    break;
     1858
     1859                  default:
     1860                    ExitFlag = 255;
     1861                    eLog() << Verbose(0) << "Invalid type given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl;
     1862                    performCriticalExit();
     1863                    break;
    16981864                }
    1699                 bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
    1700                 counter = 0;
    1701                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    1702                   Actives[counter++] = (*BigFinder)->ActiveFlag;
    1703                   (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
    1704                 }
    1705                 LCList = new LinkedCell(Boundary, 2.*radius);
    1706                 element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
    1707                 FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
    1708                 int ranges[NDIM] = {1,1,1};
    1709                 CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
    1710                 OutputCorrelationToSurface(&output, surfacemap);
    1711                 BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
    1712                 OutputCorrelation ( &binoutput, binmap );
    1713                 output.close();
    1714                 binoutput.close();
    1715                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
    1716                   (*BigFinder)->ActiveFlag = Actives[counter++];
    1717                 Free(&Actives);
    1718                 delete(LCList);
    1719                 delete(TesselStruct);
    1720                 argptr+=3;
    17211865              }
    17221866              break;
    … …  
    17371881            case 'F':
    17381882              if (ExitFlag == 0) ExitFlag = 1;
    1739               if (argptr+6 >=argc) {
     1883              MaxDistance = -1;
     1884              if (argv[argptr-1][2] == 'F') {
     1885                // fetch first argument as max distance to surface
     1886                MaxDistance = atof(argv[argptr++]);
     1887                Log() << Verbose(0) << "Filling with maximum layer distance of " << MaxDistance << "." << endl;
     1888              }
     1889              if ((argptr+7 >=argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+7]))) {
    17401890                ExitFlag = 255;
    1741                 eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
     1891                eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <xyz of filler> <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
    17421892                performCriticalExit();
    17431893              } else {
    … …  
    17461896                // construct water molecule
    17471897                molecule *filler = new molecule(periode);
     1898                if (!filler->AddXYZFile(argv[argptr])) {
     1899                  eLog() << Verbose(0) << "Could not parse filler molecule from " << argv[argptr] << "." << endl;
     1900                }
     1901                filler->SetNameFromFilename(argv[argptr]);
     1902                configuration.BG->ConstructBondGraph(filler);
    17481903                molecule *Filling = NULL;
    1749                 atom *second = NULL, *third = NULL;
    1750 //                first = new atom();
    1751 //                first->type = periode->FindElement(5);
    1752 //                first->x.Zero();
    1753 //                filler->AddAtom(first);
    1754                 first = new atom();
    1755                 first->type = periode->FindElement(1);
    1756                 first->x.Init(0.441, -0.143, 0.);
    1757                 filler->AddAtom(first);
    1758                 second = new atom();
    1759                 second->type = periode->FindElement(1);
    1760                 second->x.Init(-0.464, 1.137, 0.0);
    1761                 filler->AddAtom(second);
    1762                 third = new atom();
    1763                 third->type = periode->FindElement(8);
    1764                 third->x.Init(-0.464, 0.177, 0.);
    1765                 filler->AddAtom(third);
    1766                 filler->AddBond(first, third, 1);
    1767                 filler->AddBond(second, third, 1);
    17681904                // call routine
    17691905                double distance[NDIM];
    17701906                for (int i=0;i<NDIM;i++)
    1771                   distance[i] = atof(argv[argptr+i]);
    1772                 Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
     1907                  distance[i] = atof(argv[argptr+i+1]);
     1908                Filling = FillBoxWithMolecule(molecules, filler, configuration, MaxDistance, distance, atof(argv[argptr+4]), atof(argv[argptr+5]), atof(argv[argptr+6]), atoi(argv[argptr+7]));
    17731909                if (Filling != NULL) {
    17741910                  Filling->ActiveFlag = false;
    … …  
    17931929              }
    17941930              break;
     1931
     1932            case 'J':
     1933              if (ExitFlag == 0) ExitFlag = 1;
     1934              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1935                ExitFlag =255;
     1936                eLog() << Verbose(0) << "Missing path of adjacency file: -j <path>" << endl;
     1937                performCriticalExit();
     1938              } else {
     1939                Log() << Verbose(0) << "Storing adjacency to path " << argv[argptr] << "." << endl;
     1940                configuration.BG->ConstructBondGraph(mol);
     1941                mol->StoreAdjacencyToFile(argv[argptr]);
     1942                argptr+=1;
     1943              }
     1944              break;
     1945
     1946            case 'j':
     1947              if (ExitFlag == 0) ExitFlag = 1;
     1948              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1949                ExitFlag =255;
     1950                eLog() << Verbose(0) << "Missing path of bonds file: -j <path>" << endl;
     1951                performCriticalExit();
     1952              } else {
     1953                Log() << Verbose(0) << "Storing bonds to path " << argv[argptr] << "." << endl;
     1954                configuration.BG->ConstructBondGraph(mol);
     1955                mol->StoreBondsToFile(argv[argptr]);
     1956                argptr+=1;
     1957              }
     1958              break;
     1959
    17951960            case 'N':
    17961961              if (ExitFlag == 0) ExitFlag = 1;
    … …  
    19542119                factor[2] = atof(argv[argptr+2]);
    19552120                mol->Scale((const double ** const)&factor);
     2121                double * const cell_size = World::get()->cell_size;
    19562122                for (int i=0;i<NDIM;i++) {
    19572123                  j += i+1;
    19582124                  x.x[i] = atof(argv[NDIM+i]);
    1959                   mol->cell_size[j]*=factor[i];
     2125                  cell_size[j]*=factor[i];
    19602126                }
    19612127                delete[](factor);
    … …  
    19732139                j = -1;
    19742140                Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
     2141                double * const cell_size = World::get()->cell_size;
    19752142                for (int i=0;i<6;i++) {
    1976                   mol->cell_size[i] = atof(argv[argptr+i]);
     2143                  cell_size[i] = atof(argv[argptr+i]);
    19772144                }
    19782145                // center
    … …  
    19912158                j = -1;
    19922159                Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
     2160                double * const cell_size = World::get()->cell_size;
    19932161                for (int i=0;i<6;i++) {
    1994                   mol->cell_size[i] = atof(argv[argptr+i]);
     2162                  cell_size[i] = atof(argv[argptr+i]);
    19952163                }
    19962164                // center
    … …  
    20142182                mol->SetBoxDimension(&x);
    20152183                // translate each coordinate by boundary
     2184                double * const cell_size = World::get()->cell_size;
    20162185                j=-1;
    20172186                for (int i=0;i<NDIM;i++) {
    20182187                  j += i+1;
    20192188                  x.x[i] = atof(argv[argptr+i]);
    2020                   mol->cell_size[j] += x.x[i]*2.;
     2189                  cell_size[j] += x.x[i]*2.;
    20212190                }
    20222191                mol->Translate((const Vector *)&x);
    … …  
    21492318              } else {
    21502319                SaveFlag = true;
     2320                double * const cell_size = World::get()->cell_size;
    21512321                for (int axis = 1; axis <= NDIM; axis++) {
    21522322                  int faktor = atoi(argv[argptr++]);
    … …  
    21752345                    x.Zero();
    21762346                    y.Zero();
    2177                     y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
     2347                    y.x[abs(axis)-1] = cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
    21782348                    for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
    21792349                      x.AddVector(&y); // per factor one cell width further
    … …  
    21962366                      mol->Translate(&x);
    21972367                    }
    2198                     mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     2368                    cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    21992369                  }
    22002370                }
    … …  
    22652435  if (molecules->ListOfMolecules.size() == 0) {
    22662436    mol = new molecule(periode);
    2267     if (mol->cell_size[0] == 0.) {
     2437    double * const cell_size = World::get()->cell_size;
     2438    if (cell_size[0] == 0.) {
    22682439      Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
    22692440      for (int i=0;i<6;i++) {
    22702441        Log() << Verbose(1) << "Cell size" << i << ": ";
    2271         cin >> mol->cell_size[i];
     2442        cin >> cell_size[i];
    22722443      }
    22732444    }

  • src/config.cpp

    r85da4e rf42e95  
    1919#include "molecule.hpp"
    2020#include "periodentafel.hpp"
     21#include "World.hpp"
    2122
    2223/******************************** Functions for class ConfigFileBuffer **********************/
    … …  
    499500//        case 'j': // BoxLength
    500501//          Log() << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
     502//          double * const cell_size = World::get()->cell_size;
    501503//          for (int i=0;i<6;i++) {
    502504//            Log() << Verbose(0) << "Cell size" << i << ": ";
    503 //            cin >> mol->cell_size[i];
     505//            cin >> cell_size[i];
    504506//          }
    505507//          break;
    … …  
    689691  ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(MaxTypes), 1, critical);
    690692  if (MaxTypes == 0) {
    691     eLog() << Verbose(0) << "There are no atoms according to MaxTypes in this config file." << endl;
    692     performCriticalExit();
     693    eLog() << Verbose(1) << "There are no atoms according to MaxTypes in this config file." << endl;
     694    //performCriticalExit();
    693695  } else {
    694696    // prescan number of ions per type
    … …  
    965967  // Unit cell and magnetic field
    966968  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    967   mol->cell_size[0] = BoxLength[0];
    968   mol->cell_size[1] = BoxLength[3];
    969   mol->cell_size[2] = BoxLength[4];
    970   mol->cell_size[3] = BoxLength[6];
    971   mol->cell_size[4] = BoxLength[7];
    972   mol->cell_size[5] = BoxLength[8];
     969  double * const cell_size = World::get()->cell_size;
     970  cell_size[0] = BoxLength[0];
     971  cell_size[1] = BoxLength[3];
     972  cell_size[2] = BoxLength[4];
     973  cell_size[3] = BoxLength[6];
     974  cell_size[4] = BoxLength[7];
     975  cell_size[5] = BoxLength[8];
    973976  //if (1) fprintf(stderr,"\n");
    974977
    … …  
    11691172
    11701173  ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    1171   mol->cell_size[0] = BoxLength[0];
    1172   mol->cell_size[1] = BoxLength[3];
    1173   mol->cell_size[2] = BoxLength[4];
    1174   mol->cell_size[3] = BoxLength[6];
    1175   mol->cell_size[4] = BoxLength[7];
    1176   mol->cell_size[5] = BoxLength[8];
     1174  double * const cell_size = World::get()->cell_size;
     1175  cell_size[0] = BoxLength[0];
     1176  cell_size[1] = BoxLength[3];
     1177  cell_size[2] = BoxLength[4];
     1178  cell_size[3] = BoxLength[6];
     1179  cell_size[4] = BoxLength[7];
     1180  cell_size[5] = BoxLength[8];
    11771181  if (1) fprintf(stderr,"\n");
    11781182  config::DoPerturbation = 0;
    … …  
    13121316  // bring MaxTypes up to date
    13131317  mol->CountElements();
     1318  const double * const cell_size = World::get()->cell_size;
    13141319  ofstream * const output = new ofstream(filename, ios::out);
    13151320  if (output != NULL) {
    … …  
    13821387    *output << endl;
    13831388    *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
    1384     *output << mol->cell_size[0] << "\t" << endl;
    1385     *output << mol->cell_size[1] << "\t" << mol->cell_size[2] << "\t" << endl;
    1386     *output << mol->cell_size[3] << "\t" << mol->cell_size[4] << "\t" << mol->cell_size[5] << "\t" << endl;
     1389    *output << cell_size[0] << "\t" << endl;
     1390    *output << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
     1391    *output << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
    13871392    // FIXME
    13881393    *output << endl;

  • src/gslvector.cpp

    r85da4e rf42e95  
    66 */
    77
     8#include <cassert>
     9#include <cmath>
     10
    811#include "gslvector.hpp"
     12#include "defs.hpp"
    913
    1014/** Constructor of class GSLVector.
    … …  
    2731};
    2832
     33/** Copy constructor of class GSLVector.
     34 * Allocates GSL structures and copies components from \a *src.
     35 * \param *src source vector
     36 */
     37GSLVector::GSLVector(const GSLVector & src) : dimension(src.dimension)
     38{
     39  vector = gsl_vector_alloc(dimension);
     40  gsl_vector_memcpy (vector, src.vector);
     41};
     42
    2943/** Destructor of class GSLVector.
    3044 * Frees GSL structures
    … …  
    3347{
    3448  gsl_vector_free(vector);
    35   dimension = 0;
    3649};
    3750
    … …  
    5265 * \return  m-th element of vector
    5366 */
    54 double GSLVector::Get(size_t m)
     67double GSLVector::Get(size_t m) const
    5568{
    5669  return gsl_vector_get (vector, m);
    … …  
    7285 * \return pointer to \a m-th element
    7386 */
    74 double *GSLVector::Pointer(size_t m)
     87double *GSLVector::Pointer(size_t m) const
    7588{
    7689  return gsl_vector_ptr (vector, m);
    … …  
    8295 * \return const pointer to \a m-th element
    8396 */
    84 const double *GSLVector::const_Pointer(size_t m)
     97const double *GSLVector::const_Pointer(size_t m) const
    8598{
    8699  return gsl_vector_const_ptr (vector, m);
     100};
     101
     102/** Returns the dimension of the vector.
     103 * \return dimension of vector
     104 */
     105#ifdef HAVE_INLINE
     106inline
     107#endif
     108size_t GSLVector::GetDimension() const
     109{
     110  return dimension;
    87111};
    88112
    … …  
    128152  return gsl_vector_reverse (vector);
    129153};
     154
     155
     156/* ========================== Operators =============================== */
     157/** Compares GSLVector \a to GSLVector \a b component-wise.
     158 * \param a base GSLVector
     159 * \param b GSLVector components to add
     160 * \return a == b
     161 */
     162bool operator==(const GSLVector& a, const GSLVector& b)
     163{
     164  bool status = true;
     165  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
     166  for (size_t i=0;i<a.GetDimension();i++)
     167    status = status && (fabs(a.Get(i) - b.Get(i)) < MYEPSILON);
     168  return status;
     169};
     170
     171/** Sums GSLVector \a to this lhs component-wise.
     172 * \param a base GSLVector
     173 * \param b GSLVector components to add
     174 * \return lhs + a
     175 */
     176const GSLVector& operator+=(GSLVector& a, const GSLVector& b)
     177{
     178  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
     179  for (size_t i=0;i<a.GetDimension();i++)
     180    a.Set(i,a.Get(i)+b.Get(i));
     181  return a;
     182};
     183
     184/** Subtracts GSLVector \a from this lhs component-wise.
     185 * \param a base GSLVector
     186 * \param b GSLVector components to add
     187 * \return lhs - a
     188 */
     189const GSLVector& operator-=(GSLVector& a, const GSLVector& b)
     190{
     191  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
     192  for (size_t i=0;i<a.GetDimension();i++)
     193    a.Set(i,a.Get(i)-b.Get(i));
     194  return a;
     195};
     196
     197/** factor each component of \a a times a double \a m.
     198 * \param a base GSLVector
     199 * \param m factor
     200 * \return lhs.Get(i) * m
     201 */
     202const GSLVector& operator*=(GSLVector& a, const double m)
     203{
     204  for (size_t i=0;i<a.GetDimension();i++)
     205    a.Set(i,a.Get(i)*m);
     206  return a;
     207};
     208
     209/** Sums two GSLVectors \a  and \b component-wise.
     210 * \param a first GSLVector
     211 * \param b second GSLVector
     212 * \return a + b
     213 */
     214GSLVector const operator+(const GSLVector& a, const GSLVector& b)
     215{
     216  GSLVector x(a);
     217  for (size_t i=0;i<a.GetDimension();i++)
     218    x.Set(i,a.Get(i)+b.Get(i));
     219  return x;
     220};
     221
     222/** Subtracts GSLVector \a from \b component-wise.
     223 * \param a first GSLVector
     224 * \param b second GSLVector
     225 * \return a - b
     226 */
     227GSLVector const operator-(const GSLVector& a, const GSLVector& b)
     228{
     229  assert(a.GetDimension() == b.GetDimension() && "Dimenions of GSLVectors to compare differ");
     230  GSLVector x(a);
     231  for (size_t i=0;i<a.GetDimension();i++)
     232    x.Set(i,a.Get(i)-b.Get(i));
     233  return x;
     234};
     235
     236/** Factors given GSLVector \a a times \a m.
     237 * \param a GSLVector
     238 * \param m factor
     239 * \return m * a
     240 */
     241GSLVector const operator*(const GSLVector& a, const double m)
     242{
     243  GSLVector x(a);
     244  for (size_t i=0;i<a.GetDimension();i++)
     245    x.Set(i,a.Get(i)*m);
     246  return x;
     247};
     248
     249/** Factors given GSLVector \a a times \a m.
     250 * \param m factor
     251 * \param a GSLVector
     252 * \return m * a
     253 */
     254GSLVector const operator*(const double m, const GSLVector& a )
     255{
     256  GSLVector x(a);
     257  for (size_t i=0;i<a.GetDimension();i++)
     258    x.Set(i,a.Get(i)*m);
     259  return x;
     260};
     261
     262ostream& operator<<(ostream& ost, const GSLVector& m)
     263{
     264  ost << "(";
     265  for (size_t i=0;i<m.GetDimension();i++) {
     266    ost << m.Get(i);
     267    if (i != 2)
     268      ost << ",";
     269  }
     270  ost << ")";
     271  return ost;
     272};
     273
     274/* ====================== Checking State ============================ */
     275/** Checks whether vector has all components zero.
     276 * @return true - vector is zero, false - vector is not
     277 */
     278bool GSLVector::IsZero() const
     279{
     280  return (fabs(Get(0))+fabs(Get(1))+fabs(Get(2)) < MYEPSILON);
     281};
     282
     283/** Checks whether vector has length of 1.
     284 * @return true - vector is normalized, false - vector is not
     285 */
     286bool GSLVector::IsOne() const
     287{
     288  double NormValue = 0.;
     289  for (size_t i=dimension;--i;)
     290    NormValue += Get(i)*Get(i);
     291  return (fabs(NormValue - 1.) < MYEPSILON);
     292};
     293

  • src/gslvector.hpp

    r85da4e rf42e95  
    1818#endif
    1919
     20#include <iostream>
    2021#include <gsl/gsl_vector.h>
    2122
    … …  
    3233  GSLVector(size_t m);
    3334  GSLVector(const GSLVector * const src);
     35  GSLVector(const GSLVector & src);
    3436  ~GSLVector();
    3537
    3638  // Accessing
    3739  void SetFromDoubleArray(double *x);
    38   double Get(size_t m);
     40  double Get(size_t m) const;
    3941  void Set(size_t m, double x);
    40   double *Pointer(size_t m);
    41   const double *const_Pointer(size_t m);
     42  double *Pointer(size_t m) const;
     43  const double *const_Pointer(size_t m) const;
     44  size_t GetDimension() const;
    4245
    4346  // Initializing
    … …  
    5053  int Reverse();
    5154
     55  // checking state
     56  bool IsZero() const;
     57  bool IsOne() const;
     58
    5259private:
    5360  gsl_vector *vector;
    5461
    55   size_t dimension;
     62  const size_t dimension;
    5663};
     64
     65ostream & operator << (ostream& ost, const GSLVector &m);
     66bool operator==(const GSLVector& a, const GSLVector& b);
     67const GSLVector& operator+=(GSLVector& a, const GSLVector& b);
     68const GSLVector& operator-=(GSLVector& a, const GSLVector& b);
     69const GSLVector& operator*=(GSLVector& a, const double m);
     70GSLVector const operator*(const GSLVector& a, const double m);
     71GSLVector const operator*(const double m, const GSLVector& a);
     72GSLVector const operator+(const GSLVector& a, const GSLVector& b);
     73GSLVector const operator-(const GSLVector& a, const GSLVector& b);
     74
    5775
    5876

  • src/linkedcell.cpp

    r85da4e rf42e95  
    4646  min.Zero();
    4747  Log() << Verbose(1) << "Begin of LinkedCell" << endl;
    48   if (set->IsEmpty()) {
    49     eLog() << Verbose(1) << "set contains no linked cell nodes!" << endl;
     48  if ((set == NULL) || (set->IsEmpty())) {
     49    eLog() << Verbose(1) << "set is NULL or contains no linked cell nodes!" << endl;
    5050    return;
    5151  }

  • src/molecule.cpp

    r85da4e rf42e95  
    2323#include "tesselation.hpp"
    2424#include "vector.hpp"
     25#include "World.hpp"
    2526
    2627/************************************* Functions for class molecule *********************************/
    … …  
    4546  for(int i=MAX_ELEMENTS;i--;)
    4647    ElementsInMolecule[i] = 0;
    47   cell_size[0] = cell_size[2] = cell_size[5]= 20.;
    48   cell_size[1] = cell_size[3] = cell_size[4]= 0.;
    49   strcpy(name,"none");
     48  strcpy(name,World::get()->DefaultName);
    5049};
    5150
    … …  
    159158  double *matrix = NULL;
    160159  bond *Binder = NULL;
     160  double * const cell_size = World::get()->cell_size;
    161161
    162162//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
    … …  
    603603void molecule::SetBoxDimension(Vector *dim)
    604604{
     605  double * const cell_size = World::get()->cell_size;
    605606  cell_size[0] = dim->x[0];
    606607  cell_size[1] = 0.;
    … …  
    693694bool molecule::CheckBounds(const Vector *x) const
    694695{
     696  double * const cell_size = World::get()->cell_size;
    695697  bool result = true;
    696698  int j =-1;

  • src/molecule.hpp

    r85da4e rf42e95  
    8282class molecule : public PointCloud {
    8383  public:
    84     double cell_size[6];//!< cell size
    8584    const periodentafel * const elemente; //!< periodic table with each element
    8685    atom *start;        //!< start of atom list
    … …  
    269268  int FragmentMolecule(int Order, config *configuration);
    270269  bool CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
     270  bool StoreBondsToFile(char *path);
    271271  bool StoreAdjacencyToFile(char *path);
    272272  bool CheckAdjacencyFileAgainstMolecule(char *path, atom **ListOfAtoms);

  • src/molecule_fragmentation.cpp

    r85da4e rf42e95  
    1818#include "molecule.hpp"
    1919#include "periodentafel.hpp"
     20#include "World.hpp"
    2021
    2122/************************************* Functions for class molecule *********************************/
    … …  
    843844
    844845  Leaf->BondDistance = mol->BondDistance;
    845   for(int i=NDIM*2;i--;)
    846     Leaf->cell_size[i] = mol->cell_size[i];
    847846
    848847  // first create the minimal set of atoms from the KeySet
    … …  
    16541653  atom *Walker = NULL;
    16551654  atom *OtherWalker = NULL;
     1655  double * const cell_size = World::get()->cell_size;
    16561656  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    16571657  enum Shading *ColorList = NULL;

  • src/molecule_geometry.cpp

    r85da4e rf42e95  
    1515#include "memoryallocator.hpp"
    1616#include "molecule.hpp"
     17#include "World.hpp"
    1718
    1819/************************************* Functions for class molecule *********************************/
    … …  
    2627  bool status = true;
    2728  const Vector *Center = DetermineCenterOfAll();
     29  double * const cell_size = World::get()->cell_size;
    2830  double *M = ReturnFullMatrixforSymmetric(cell_size);
    2931  double *Minv = InverseMatrix(M);
    … …  
    4648{
    4749  bool status = true;
     50  double * const cell_size = World::get()->cell_size;
    4851  double *M = ReturnFullMatrixforSymmetric(cell_size);
    4952  double *Minv = InverseMatrix(M);
    … …  
    226229void molecule::TranslatePeriodically(const Vector *trans)
    227230{
     231  double * const cell_size = World::get()->cell_size;
    228232  double *M = ReturnFullMatrixforSymmetric(cell_size);
    229233  double *Minv = InverseMatrix(M);
    … …  
    252256{
    253257  atom *Walker = start;
     258  double * const cell_size = World::get()->cell_size;
    254259  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    255260  double *inversematrix = InverseMatrix(cell_size);

  • src/molecule_graph.cpp

    r85da4e rf42e95  
    1717#include "memoryallocator.hpp"
    1818#include "molecule.hpp"
     19#include "World.hpp"
    1920
    2021struct BFSAccounting
    … …  
    106107  LinkedCell *LC = NULL;
    107108  bool free_BG = false;
     109  double * const cell_size = World::get()->cell_size;
    108110
    109111  if (BG == NULL) {
    … …  
    493495  bond *Binder = NULL;
    494496
     497  if (AtomCount == 0)
     498    return SubGraphs;
    495499  Log() << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;
    496500  DepthFirstSearchAnalysis_Init(DFS, this);
    … …  
    992996  Log() << Verbose(1) << "Saving adjacency list ... ";
    993997  if (AdjacencyFile != NULL) {
     998    AdjacencyFile << "m\tn" << endl;
    994999    ActOnAllAtoms(&atom::OutputAdjacency, &AdjacencyFile);
    9951000    AdjacencyFile.close();
     1001    Log() << Verbose(1) << "done." << endl;
     1002  } else {
     1003    Log() << Verbose(1) << "failed to open file " << line.str() << "." << endl;
     1004    status = false;
     1005  }
     1006
     1007  return status;
     1008}
     1009;
     1010
     1011/** Storing the bond structure of a molecule to file.
     1012 * Simply stores Atom::nr and then the Atom::nr of all bond partners, one per line.
     1013 * \param *out output stream for debugging
     1014 * \param *path path to file
     1015 * \return true - file written successfully, false - writing failed
     1016 */
     1017bool molecule::StoreBondsToFile(char *path)
     1018{
     1019  ofstream BondFile;
     1020  stringstream line;
     1021  bool status = true;
     1022
     1023  line << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
     1024  BondFile.open(line.str().c_str(), ios::out);
     1025  Log() << Verbose(1) << "Saving adjacency list ... ";
     1026  if (BondFile != NULL) {
     1027    BondFile << "m\tn" << endl;
     1028    ActOnAllAtoms(&atom::OutputBonds, &BondFile);
     1029    BondFile.close();
    9961030    Log() << Verbose(1) << "done." << endl;
    9971031  } else {

  • src/moleculelist.cpp

    r85da4e rf42e95  
    1919#include "memoryallocator.hpp"
    2020#include "periodentafel.hpp"
     21#include "World.hpp"
    2122
    2223/*********************************** Functions for class MoleculeListClass *************************/
    … …  
    638639  int FragmentCounter = 0;
    639640  ofstream output;
    640 
     641  double cell_size_backup[6];
     642  double * const cell_size = World::get()->cell_size;
     643
     644  // backup cell_size
     645  for (int i=0;i<6;i++)
     646    cell_size_backup[i] = cell_size[i];
    641647  // store the fragments as config and as xyz
    642648  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    … …  
    682688      j += k + 1;
    683689      BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
    684       (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
     690      cell_size[j] = BoxDimension.x[k] * 2.;
    685691    }
    686692    (*ListRunner)->Translate(&BoxDimension);
    … …  
    724730  // printing final number
    725731  Log() << Verbose(2) << "Final number of fragments: " << FragmentCounter << "." << endl;
     732
     733  // restore cell_size
     734  for (int i=0;i<6;i++)
     735    cell_size[i] = cell_size_backup[i];
    726736
    727737  return result;
    … …  
    776786
    777787  // 1. dissect the molecule into connected subgraphs
    778   configuration->BG->ConstructBondGraph(mol);
     788  if (configuration->BG->ConstructBondGraph(mol)) {
     789    delete (mol);
     790    eLog() << Verbose(1) << "There are no bonds." << endl;
     791    return;
     792  }
    779793
    780794  // 2. scan for connected subgraphs
    … …  
    783797  Subgraphs = mol->DepthFirstSearchAnalysis(BackEdgeStack);
    784798  delete(BackEdgeStack);
     799  if ((Subgraphs == NULL) || (Subgraphs->next == NULL)) {
     800    delete (mol);
     801    eLog() << Verbose(1) << "There are no atoms." << endl;
     802    return;
     803  }
    785804
    786805  // 3. dissect (the following construct is needed to have the atoms not in the order of the DFS, but in

  • src/tesselation.cpp

    r85da4e rf42e95  
    1414#include "tesselation.hpp"
    1515#include "tesselationhelpers.hpp"
     16#include "triangleintersectionlist.hpp"
    1617#include "vector.hpp"
    1718#include "verbose.hpp"
    … …  
    474475};
    475476
    476 /** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses.
    477  * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane
     477/** Finds the point on the triangle to the point \a *x.
     478 * We call Vector::GetIntersectionWithPlane() with \a * and the center of the triangle to receive an intersection point.
     479 * Then we check the in-plane part (the part projected down onto plane). We check whether it crosses one of the
     480 * boundary lines. If it does, we return this intersection as closest point, otherwise the projected point down.
    478481 * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not.
    479482 * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line
    … …  
    32223225};
    32233226
    3224 
    32253227/** Finds the triangle that is closest to a given Vector \a *x.
    32263228 * \param *out output stream for debugging
    … …  
    33153317 * \param *out output stream for debugging
    33163318 * \param *x Vector to look from
     3319 * \param &distance contains found distance on return
    33173320 * \return list of BoundaryTriangleSet of nearest triangles or NULL.
    33183321 */
    … …  
    33583361bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
    33593362{
    3360   return (GetDistanceSquaredToSurface(Point, LC) < MYEPSILON);
    3361 }
     3363  Info FunctionInfo(__func__);
     3364  TriangleIntersectionList Intersections(&Point,this,LC);
     3365
     3366  return Intersections.IsInside();
     3367};
    33623368
    33633369/** Returns the distance to the surface given by the tesselation.
    … …  
    34323438};
    34333439
    3434 /** Calculates distance to a tesselated surface.
     3440/** Calculates minimum distance from \a&Point to a tesselated surface.
    34353441 * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
    34363442 * \param &Point point to calculate distance from
    … …  
    34383444 * \return distance squared to closest point on surface
    34393445 */
    3440 double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const
    3441 {
    3442   BoundaryTriangleSet *triangle = FindClosestTriangleToVector(&Point, LC);
    3443   const double distance = GetDistanceSquaredToTriangle(Point, triangle);
    3444   return Min(distance, LC->RADIUS);
     3446double Tesselation::GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const
     3447{
     3448  Info FunctionInfo(__func__);
     3449  TriangleIntersectionList Intersections(&Point,this,LC);
     3450
     3451  return Intersections.GetSmallestDistance();
     3452};
     3453
     3454/** Calculates minimum distance from \a&Point to a tesselated surface.
     3455 * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
     3456 * \param &Point point to calculate distance from
     3457 * \param *LC needed for finding closest points fast
     3458 * \return distance squared to closest point on surface
     3459 */
     3460BoundaryTriangleSet * Tesselation::GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const
     3461{
     3462  Info FunctionInfo(__func__);
     3463  TriangleIntersectionList Intersections(&Point,this,LC);
     3464
     3465  return Intersections.GetClosestTriangle();
    34453466};
    34463467

  • src/tesselation.hpp

    r85da4e rf42e95  
    318318    bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const;
    319319    double GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) const;
    320     double GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const;
     320    double GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const;
     321    BoundaryTriangleSet * GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const;
    321322    bool AddBoundaryPoint(TesselPoint * Walker, const int n);
    322323    DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;

  • src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

    r85da4e rf42e95  
    6060
    6161  // construct molecule (tetraeder of hydrogens) base
     62  TestSurfaceMolecule = new molecule(tafel);
     63  Walker = new atom();
     64  Walker->type = hydrogen;
     65  Walker->node->Init(1., 0., 1. );
     66  TestSurfaceMolecule->AddAtom(Walker);
     67  Walker = new atom();
     68  Walker->type = hydrogen;
     69  Walker->node->Init(0., 1., 1. );
     70  TestSurfaceMolecule->AddAtom(Walker);
     71  Walker = new atom();
     72  Walker->type = hydrogen;
     73  Walker->node->Init(1., 1., 0. );
     74  TestSurfaceMolecule->AddAtom(Walker);
     75  Walker = new atom();
     76  Walker->type = hydrogen;
     77  Walker->node->Init(0., 0., 0. );
     78  TestSurfaceMolecule->AddAtom(Walker);
     79
     80  // check that TestMolecule was correctly constructed
     81  CPPUNIT_ASSERT_EQUAL( TestSurfaceMolecule->AtomCount, 4 );
     82
     83  TestList = new MoleculeListClass;
     84  TestSurfaceMolecule->ActiveFlag = true;
     85  TestList->insert(TestSurfaceMolecule);
     86
     87  // init tesselation and linked cell
     88  Surface = new Tesselation;
     89  LC = new LinkedCell(TestSurfaceMolecule, 5.);
     90  FindNonConvexBorder(TestSurfaceMolecule, Surface, (const LinkedCell *&)LC, 2.5, NULL);
     91
     92  // add outer atoms
    6293  TestMolecule = new molecule(tafel);
    6394  Walker = new atom();
    64   Walker->type = hydrogen;
    65   Walker->node->Init(1., 0., 1. );
    66   TestMolecule->AddAtom(Walker);
    67   Walker = new atom();
    68   Walker->type = hydrogen;
    69   Walker->node->Init(0., 1., 1. );
    70   TestMolecule->AddAtom(Walker);
    71   Walker = new atom();
    72   Walker->type = hydrogen;
    73   Walker->node->Init(1., 1., 0. );
    74   TestMolecule->AddAtom(Walker);
    75   Walker = new atom();
    76   Walker->type = hydrogen;
    77   Walker->node->Init(0., 0., 0. );
    78   TestMolecule->AddAtom(Walker);
    79 
    80   // check that TestMolecule was correctly constructed
    81   CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 );
    82 
    83   TestList = new MoleculeListClass;
     95  Walker->type = carbon;
     96  Walker->node->Init(4., 0., 4. );
     97  TestMolecule->AddAtom(Walker);
     98  Walker = new atom();
     99  Walker->type = carbon;
     100  Walker->node->Init(0., 4., 4. );
     101  TestMolecule->AddAtom(Walker);
     102  Walker = new atom();
     103  Walker->type = carbon;
     104  Walker->node->Init(4., 4., 0. );
     105  TestMolecule->AddAtom(Walker);
     106  // add inner atoms
     107  Walker = new atom();
     108  Walker->type = carbon;
     109  Walker->node->Init(0.5, 0.5, 0.5 );
     110  TestMolecule->AddAtom(Walker);
    84111  TestMolecule->ActiveFlag = true;
    85112  TestList->insert(TestMolecule);
    86 
    87   // init tesselation and linked cell
    88   Surface = new Tesselation;
    89   FindNonConvexBorder(TestMolecule, Surface, (const LinkedCell *&)LC, 2.5, NULL);
    90   LC = new LinkedCell(TestMolecule, 5.);
    91   CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->PointsOnBoundary.size() );
    92   CPPUNIT_ASSERT_EQUAL( (size_t)6, Surface->LinesOnBoundary.size() );
    93   CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->TrianglesOnBoundary.size() );
    94 
    95   // add outer atoms
    96   Walker = new atom();
    97   Walker->type = carbon;
    98   Walker->node->Init(4., 0., 4. );
    99   TestMolecule->AddAtom(Walker);
    100   Walker = new atom();
    101   Walker->type = carbon;
    102   Walker->node->Init(0., 4., 4. );
    103   TestMolecule->AddAtom(Walker);
    104   Walker = new atom();
    105   Walker->type = carbon;
    106   Walker->node->Init(4., 4., 0. );
    107   TestMolecule->AddAtom(Walker);
    108   // add inner atoms
    109   Walker = new atom();
    110   Walker->type = carbon;
    111   Walker->node->Init(0.5, 0.5, 0.5 );
    112   TestMolecule->AddAtom(Walker);
    113113
    114114  // init maps
    … …  
    136136
    137137
     138/** Checks whether setup() does the right thing.
     139 */
     140void AnalysisCorrelationToSurfaceUnitTest::SurfaceTest()
     141{
     142  CPPUNIT_ASSERT_EQUAL( 4, TestSurfaceMolecule->AtomCount );
     143  CPPUNIT_ASSERT_EQUAL( 4, TestMolecule->AtomCount );
     144  CPPUNIT_ASSERT_EQUAL( (size_t)2, TestList->ListOfMolecules.size() );
     145  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->PointsOnBoundary.size() );
     146  CPPUNIT_ASSERT_EQUAL( (size_t)6, Surface->LinesOnBoundary.size() );
     147  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->TrianglesOnBoundary.size() );
     148};
     149
    138150void AnalysisCorrelationToSurfaceUnitTest::CorrelationToSurfaceTest()
    139151{
    140152  // do the pair correlation
    141153  surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
     154//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    142155  CPPUNIT_ASSERT( surfacemap != NULL );
    143156  CPPUNIT_ASSERT_EQUAL( (size_t)4, surfacemap->size() );
    … …  
    149162  surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
    150163  // put pair correlation into bins and check with no range
     164//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    151165  binmap = BinData( surfacemap, 0.5, 0., 0. );
    152166  CPPUNIT_ASSERT_EQUAL( (size_t)1, binmap->size() );
    153   //OutputCorrelation ( binmap );
     167  OutputCorrelation ( (ofstream *)&cout, binmap );
    154168  tester = binmap->begin();
    155169  CPPUNIT_ASSERT_EQUAL( 0., tester->first );
    … …  
    162176  BinPairMap::iterator tester;
    163177  surfacemap = CorrelationToSurface( TestList, hydrogen, Surface, LC );
     178//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    164179  // ... and check with [0., 2.] range
    165180  binmap = BinData( surfacemap, 0.5, 0., 2. );
    166181  CPPUNIT_ASSERT_EQUAL( (size_t)5, binmap->size() );
    167   //OutputCorrelation ( binmap );
     182//  OutputCorrelation ( (ofstream *)&cout, binmap );
    168183  tester = binmap->begin();
    169184  CPPUNIT_ASSERT_EQUAL( 0., tester->first );
    … …  
    179194  BinPairMap::iterator tester;
    180195  surfacemap = CorrelationToSurface( TestList, carbon, Surface, LC );
     196//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    181197  // put pair correlation into bins and check with no range
    182198  binmap = BinData( surfacemap, 0.5, 0., 0. );
    … …  
    184200  OutputCorrelation ( (ofstream *)&cout, binmap );
    185201  // inside point is first and must have negative value
    186   tester = binmap->lower_bound(2.95); // start depends on the min value and
     202  tester = binmap->lower_bound(4.25-0.5); // start depends on the min value and
    187203  CPPUNIT_ASSERT( tester != binmap->end() );
    188204  CPPUNIT_ASSERT_EQUAL( 3, tester->second );
    189205  // inner point
    190   tester = binmap->lower_bound(-0.5);
     206  tester = binmap->lower_bound(0.);
    191207  CPPUNIT_ASSERT( tester != binmap->end() );
    192208  CPPUNIT_ASSERT_EQUAL( 1, tester->second );
    … …  
    197213  BinPairMap::iterator tester;
    198214  surfacemap = CorrelationToSurface( TestList, carbon, Surface, LC );
     215//  OutputCorrelationToSurface ( (ofstream *)&cout, surfacemap );
    199216  // ... and check with [0., 2.] range
    200217  binmap = BinData( surfacemap, 0.5, -2., 4. );
    201218  CPPUNIT_ASSERT_EQUAL( (size_t)13, binmap->size() );
    202   OutputCorrelation ( (ofstream *)&cout, binmap );
     219//  OutputCorrelation ( (ofstream *)&cout, binmap );
    203220  // three outside points
    204   tester = binmap->lower_bound(3.);
     221  tester = binmap->lower_bound(4.25-0.5);
    205222  CPPUNIT_ASSERT( tester != binmap->end() );
    206223  CPPUNIT_ASSERT_EQUAL( 3, tester->second );
    207224  // inner point
    208   tester = binmap->lower_bound(-0.5);
     225  tester = binmap->lower_bound(0.);
    209226  CPPUNIT_ASSERT( tester != binmap->end() );
    210227  CPPUNIT_ASSERT_EQUAL( 1, tester->second );
    211 
    212228};
    213229

  • src/unittests/AnalysisCorrelationToSurfaceUnitTest.hpp

    r85da4e rf42e95  
    2323{
    2424    CPPUNIT_TEST_SUITE( AnalysisCorrelationToSurfaceUnitTest ) ;
     25    CPPUNIT_TEST ( SurfaceTest );
    2526    CPPUNIT_TEST ( CorrelationToSurfaceTest );
    2627    CPPUNIT_TEST ( CorrelationToSurfaceHydrogenBinNoRangeTest );
    … …  
    3334      void setUp();
    3435      void tearDown();
     36      void SurfaceTest();
    3537      void CorrelationToSurfaceTest();
    3638      void CorrelationToSurfaceHydrogenBinNoRangeTest();
    … …  
    4345      MoleculeListClass *TestList;
    4446      molecule *TestMolecule;
     47      molecule *TestSurfaceMolecule;
    4548      element *hydrogen;
    4649      element *carbon;

  • src/unittests/gslvectorunittest.cpp

    r85da4e rf42e95  
    115115
    116116
     117/** UnitTest for operators.
     118 */
     119void GSLVectorTest::OperatorIsTest()
     120{
     121  GSLVector zero(3);
     122  GSLVector unit(3);
     123  zero.SetZero();
     124  unit.SetZero();
     125  unit.Set(1,1.);
     126  // summation and scaling
     127  CPPUNIT_ASSERT_EQUAL( true, unit.IsOne() );
     128  CPPUNIT_ASSERT_EQUAL( false, zero.IsOne() );
     129  CPPUNIT_ASSERT_EQUAL( false, unit.IsZero() );
     130  CPPUNIT_ASSERT_EQUAL( true, zero.IsZero() );
     131};
     132
     133/** UnitTest for operators.
     134 */
     135void GSLVectorTest::OperatorAlgebraTest()
     136{
     137  GSLVector zero(3);
     138  GSLVector unit(3);
     139  zero.SetZero();
     140  unit.SetZero();
     141  unit.Set(1,1.);
     142  // summation and scaling
     143  CPPUNIT_ASSERT_EQUAL( true, (zero+unit).IsOne() );
     144  CPPUNIT_ASSERT_EQUAL( true, (zero+unit).IsOne() );
     145  CPPUNIT_ASSERT_EQUAL( true, (zero-unit).IsOne() );
     146  CPPUNIT_ASSERT_EQUAL( false, (zero-unit).IsZero() );
     147  CPPUNIT_ASSERT_EQUAL( true, (zero+zero).IsZero() );
     148  CPPUNIT_ASSERT_EQUAL( false, (unit*0.98).IsOne() );
     149  CPPUNIT_ASSERT_EQUAL( false, (0.98*unit).IsOne() );
     150  CPPUNIT_ASSERT_EQUAL( true, (unit*1.).IsOne() );
     151  CPPUNIT_ASSERT_EQUAL( true, (1.*unit).IsOne() );
     152
     153  CPPUNIT_ASSERT_EQUAL( unit, (zero+unit) );
     154  CPPUNIT_ASSERT_EQUAL( zero, (zero+zero) );
     155  CPPUNIT_ASSERT_EQUAL( unit, (unit+zero) );
     156
     157  unit += zero;
     158  CPPUNIT_ASSERT_EQUAL( true, unit.IsOne() );
     159  unit *= 1.;
     160  CPPUNIT_ASSERT_EQUAL( true, unit.IsOne() );
     161};
     162
    117163/********************************************** Main routine **************************************/
    118164

  • src/unittests/gslvectorunittest.hpp

    r85da4e rf42e95  
    2222    CPPUNIT_TEST (CopyTest );
    2323    CPPUNIT_TEST (ExchangeTest );
     24    CPPUNIT_TEST (OperatorAlgebraTest );
     25    CPPUNIT_TEST (OperatorIsTest );
    2426    CPPUNIT_TEST_SUITE_END();
    2527
    … …  
    3234    void CopyTest();
    3335    void ExchangeTest();
     36    void OperatorIsTest();
     37    void OperatorAlgebraTest();
    3438
    3539private:

  • src/vector.cpp

    r85da4e rf42e95  
    1515#include "vector.hpp"
    1616#include "verbose.hpp"
     17#include "World.hpp"
    1718
    1819#include <gsl/gsl_linalg.h>
    … …  
    2627 */
    2728Vector::Vector() { x[0] = x[1] = x[2] = 0.; };
     29
     30/** Constructor of class vector.
     31 */
     32Vector::Vector(const Vector * const a)
     33{
     34  x[0] = a->x[0];
     35  x[1] = a->x[1];
     36  x[2] = a->x[2];
     37};
     38
     39/** Constructor of class vector.
     40 */
     41Vector::Vector(const Vector &a)
     42{
     43  x[0] = a.x[0];
     44  x[1] = a.x[1];
     45  x[2] = a.x[2];
     46};
    2847
    2948/** Constructor of class vector.
    … …  
    267286};
    268287
    269 /** Calculates the minimum distance of this vector to the plane.
     288/** Calculates the minimum distance vector of this vector to the plane.
    270289 * \param *out output stream for debugging
    271290 * \param *PlaneNormal normal of plane
    272291 * \param *PlaneOffset offset of plane
    273  * \return distance to plane
    274  */
    275 double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
     292 * \return distance vector onto to plane
     293 */
     294Vector Vector::GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
    276295{
    277296  Vector temp;
    … …  
    291310    sign = 0.;
    292311
    293   return (temp.Norm()*sign);
     312  temp.Normalize();
     313  temp.Scale(sign);
     314  return temp;
     315};
     316
     317/** Calculates the minimum distance of this vector to the plane.
     318 * \sa Vector::GetDistanceVectorToPlane()
     319 * \param *out output stream for debugging
     320 * \param *PlaneNormal normal of plane
     321 * \param *PlaneOffset offset of plane
     322 * \return distance to plane
     323 */
     324double Vector::DistanceToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const
     325{
     326  return GetDistanceVectorToPlane(PlaneNormal,PlaneOffset).Norm();
    294327};
    295328

  • src/vector.hpp

    r85da4e rf42e95  
    2727
    2828  Vector();
     29  Vector(const Vector * const a);
     30  Vector(const Vector &a);
    2931  Vector(const double x1, const double x2, const double x3);
    3032  ~Vector();
    … …  
    6769  void LinearCombinationOfVectors(const Vector * const x1, const Vector * const x2, const Vector * const x3, const double * const factors);
    6870  double CutsPlaneAt(const Vector * const A, const Vector * const B, const Vector * const C) const;
     71  Vector GetDistanceVectorToPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset) const;
    6972  bool GetIntersectionWithPlane(const Vector * const PlaneNormal, const Vector * const PlaneOffset, const Vector * const Origin, const Vector * const LineVector);
    7073  bool GetIntersectionOfTwoLinesOnPlane(const Vector * const Line1a, const Vector * const Line1b, const Vector * const Line2a, const Vector * const Line2b, const Vector *Normal = NULL);
Note: See TracChangeset for help on using the changeset viewer.