Changes in / [a7b761b:8f215d]


Ignore:
Location:
src
Files:
1 added
1 deleted
39 edited

Legend:

Unmodified
Added
Removed
  • src/Helpers/Assert.cpp

    ra7b761b r8f215d  
    5454
    5555
     56
    5657bool _my_assert::check(const bool res,
    5758                       const char* condition,
     
    6263{
    6364  if(!res){
    64     cout << "Assertion \"" << condition << "\" failed in file " << filename << " at line " << line << endl;
     65    cout << "Assertion " << condition << " failed in file " << filename << " at line " << line << endl;
    6566    cout << "Assertion Message: " << message << std::endl;
    6667    while(true){
  • src/Legacy/oldmenu.cpp

    ra7b761b r8f215d  
    423423void oldmenu::RemoveAtoms(molecule *mol)
    424424{
    425   atom *second;
     425  atom *first, *second;
    426426  int axis;
    427427  double tmp1, tmp2;
     
    446446      break;
    447447    case 'b':
    448       {
    449         second = mol->AskAtom("Enter number of atom as reference point: ");
    450         Log() << Verbose(0) << "Enter radius: ";
    451         cin >> tmp1;
    452         molecule::iterator runner;
    453         for (molecule::iterator iter = mol->begin(); iter != mol->end(); ) {
    454           runner = iter++;
    455           if ((*runner)->x.DistanceSquared((*runner)->x) > tmp1*tmp1) // distance to first above radius ...
    456             mol->RemoveAtom((*runner));
    457         }
     448      second = mol->AskAtom("Enter number of atom as reference point: ");
     449      Log() << Verbose(0) << "Enter radius: ";
     450      cin >> tmp1;
     451      first = mol->start;
     452      second = first->next;
     453      while(second != mol->end) {
     454        first = second;
     455        second = first->next;
     456        if (first->x.DistanceSquared(second->x) > tmp1*tmp1) // distance to first above radius ...
     457          mol->RemoveAtom(first);
    458458      }
    459459      break;
     
    465465      Log() << Verbose(0) << "Upper boundary: ";
    466466      cin >> tmp2;
    467       molecule::iterator runner;
    468       for (molecule::iterator iter = mol->begin(); iter != mol->end(); ) {
    469         runner = iter++;
    470         if (((*runner)->x[axis] < tmp1) || ((*runner)->x[axis] > tmp2)) {// out of boundary ...
    471           //Log() << Verbose(0) << "Atom " << *(*runner) << " with " << (*runner)->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
    472           mol->RemoveAtom((*runner));
     467      first = mol->start;
     468      second = first->next;
     469      while(second != mol->end) {
     470        first = second;
     471        second = first->next;
     472        if ((first->x[axis] < tmp1) || (first->x[axis] > tmp2)) {// out of boundary ...
     473          //Log() << Verbose(0) << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
     474          mol->RemoveAtom(first);
    473475        }
    474476      }
     
    513515        min[i] = 0.;
    514516
    515       for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    516         Z = (*iter)->type->Z;
     517      second = mol->start;
     518      while ((second->next != mol->end)) {
     519        second = second->next; // advance
     520        Z = second->type->Z;
    517521        tmp1 = 0.;
    518         if (first != (*iter)) {
    519           x = first->x - (*iter)->x;
     522        if (first != second) {
     523          x = first->x - second->x;
    520524          tmp1 = x.Norm();
    521525        }
    522526        if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
    523         //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << ((*iter)->nr << ": " << tmp1 << " a.u." << endl;
     527        //Log() << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
    524528      }
    525529      for (int i=MAX_ELEMENTS;i--;)
     
    750754    Log() << Verbose(0) << "State the factor: ";
    751755    cin >> faktor;
    752     if (mol->getAtomCount() != 0) {  // if there is more than none
    753       count = mol->getAtomCount();  // is changed becausing of adding, thus has to be stored away beforehand
     756
     757    mol->CountAtoms(); // recount atoms
     758    if (mol->AtomCount != 0) {  // if there is more than none
     759      count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
    754760      Elements = new const element *[count];
    755761      vectors = new Vector *[count];
    756762      j = 0;
    757       for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    758         Elements[j] = (*iter)->type;
    759         vectors[j] = &(*iter)->x;
     763      first = mol->start;
     764      while (first->next != mol->end) { // make a list of all atoms with coordinates and element
     765        first = first->next;
     766        Elements[j] = first->type;
     767        vectors[j] = &first->x;
    760768        j++;
    761769      }
     
    10161024    return;
    10171025  }
     1026  atom *Walker = mol->start;
    10181027
    10191028  // generate some KeySets
    10201029  Log() << Verbose(0) << "Generating KeySets." << endl;
    1021   KeySet TestSets[mol->getAtomCount()+1];
     1030  KeySet TestSets[mol->AtomCount+1];
    10221031  i=1;
    1023   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     1032  while (Walker->next != mol->end) {
     1033    Walker = Walker->next;
    10241034    for (int j=0;j<i;j++) {
    1025       TestSets[j].insert((*iter)->nr);
     1035      TestSets[j].insert(Walker->nr);
    10261036    }
    10271037    i++;
     
    10291039  Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl;
    10301040  KeySetTestPair test;
    1031   molecule::const_iterator iter = mol->begin();
    1032   if (iter != mol->end()) {
    1033     test = TestSets[mol->getAtomCount()-1].insert((*iter)->nr);
    1034     if (test.second) {
    1035       Log() << Verbose(1) << "Insertion worked?!" << endl;
    1036     } else {
    1037       Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
    1038     }
     1041  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
     1042  if (test.second) {
     1043    Log() << Verbose(1) << "Insertion worked?!" << endl;
    10391044  } else {
    1040     eLog() << Verbose(1) << "No atoms to test double insertion." << endl;
    1041   }
     1045    Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
     1046  }
     1047  TestSets[mol->AtomCount].insert(mol->end->previous->nr);
     1048  TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
    10421049
    10431050  // constructing Graph structure
     
    10471054  // insert KeySets into Subgraphs
    10481055  Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl;
    1049   for (int j=0;j<mol->getAtomCount();j++) {
     1056  for (int j=0;j<mol->AtomCount;j++) {
    10501057    Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
    10511058  }
    10521059  Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl;
    10531060  GraphTestPair test2;
    1054   test2 = Subgraphs.insert(GraphPair (TestSets[mol->getAtomCount()],pair<int, double>(counter++, 1.)));
     1061  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
    10551062  if (test2.second) {
    10561063    Log() << Verbose(1) << "Insertion worked?!" << endl;
  • src/Makefile.am

    ra7b761b r8f215d  
    7676                                   Descriptors/MoleculeIdDescriptor.cpp
    7777                                   
    78 
     78                                   
    7979DESCRIPTORHEADER = Descriptors/AtomDescriptor.hpp \
    8080                                   Descriptors/AtomIdDescriptor.hpp \
     
    9494                                  Exceptions/MathException.hpp \
    9595                                  Exceptions/ZeroVectorException.hpp
    96 
    9796
    9897SOURCE = ${ANALYSISSOURCE} \
     
    113112                 graph.cpp \
    114113                 helpers.cpp \
    115                  Helpers/Assert.cpp \
    116114                 info.cpp \
    117115                 leastsquaremin.cpp \
    118116                 Line.cpp \
    119117                 linkedcell.cpp \
     118                 lists.cpp \
    120119                 log.cpp \
    121120                 logger.cpp \
     
    134133                 tesselationhelpers.cpp \
    135134                 triangleintersectionlist.cpp \
    136                  vector.cpp \
    137135                 verbose.cpp \
    138136                 vector_ops.cpp \
  • src/Patterns/Cacheable.hpp

    ra7b761b r8f215d  
    2828        owner(_owner)
    2929        {}
    30       virtual T& getValue()=0;
     30      virtual T getValue()=0;
    3131      virtual void invalidate()=0;
    3232      virtual bool isValid()=0;
     
    4646        {}
    4747
    48       virtual T& getValue(){
     48      virtual T getValue(){
    4949        // set the state to valid
    5050        State::owner->switchState(State::owner->validState);
     
    7272        {}
    7373
    74       virtual T& getValue(){
     74      virtual T getValue(){
    7575        return content;
    7676      }
     
    100100        {}
    101101
    102       virtual T& getValue(){
     102      virtual T getValue(){
    103103        ASSERT(0,"Cannot get a value from a Cacheable after it's Observable has died");
    104104        // we have to return a grossly invalid reference, because no value can be produced anymore
     
    134134    void subjectKilled(Observable *subject);
    135135  private:
     136
    136137    void switchState(state_ptr newState);
    137138
     
    143144
    144145    Observable *owner;
     146
    145147    boost::function<T()> recalcMethod;
    146148
     
    219221
    220222    const bool isValid() const;
    221     const T& operator*() const;
     223    const T operator*() const;
    222224
    223225    // methods implemented for base-class Observer
     
    235237
    236238  template<typename T>
    237   const T& Cacheable<T>::operator*() const{
     239  const T Cacheable<T>::operator*() const{
    238240    return recalcMethod();
    239241  }
  • src/Patterns/Observer.cpp

    ra7b761b r8f215d  
    8282Observable::_Observable_protector::_Observable_protector(Observable *_protege) :
    8383  protege(_protege)
    84 {
    85   start_observer_internal(protege);
    86 }
    87 
    88 Observable::_Observable_protector::_Observable_protector(const _Observable_protector &dest) :
    89     protege(dest.protege)
    9084{
    9185  start_observer_internal(protege);
     
    195189  ASSERT(callTable.count(this),"SignOff called for an Observable without Observers.");
    196190  callees_t &callees = callTable[this];
    197 
    198191  callees_t::iterator iter;
    199192  callees_t::iterator deliter;
  • src/Patterns/Observer.hpp

    ra7b761b r8f215d  
    3636typedef Notification *const Notification_ptr;
    3737
    38 template<class _Set>
    39 class ObservedIterator;
    40 
    4138/**
    4239 * An Observer is notified by all Observed objects, when anything changes.
     
    5653  friend class Observable;
    5754  friend class Notification;
    58   template<class> friend class ObservedIterator;
    59 
    6055public:
    6156  Observer();
     
    157152  static std::set<Observable*> busyObservables;
    158153
     154
    159155  //! @cond
    160156  // Structure for RAII-Style notification
     
    168164  public:
    169165    _Observable_protector(Observable *);
    170     _Observable_protector(const _Observable_protector&);
    171166    ~_Observable_protector();
    172167  private:
  • src/analysis_bonds.cpp

    ra7b761b r8f215d  
    2626  Mean = 0.;
    2727
     28  atom *Walker = mol->start;
    2829  int AtomCount = 0;
    29   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    30     const int count = (*iter)->ListOfBonds.size();
     30  while (Walker->next != mol->end) {
     31    Walker = Walker->next;
     32    const int count = Walker->ListOfBonds.size();
    3133    if (Max < count)
    3234      Max = count;
     
    5658
    5759  int AtomNo = 0;
    58   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    59     if ((*iter)->type == type1)
    60       for (BondList::const_iterator BondRunner = (*iter)->ListOfBonds.begin(); BondRunner != (*iter)->ListOfBonds.end(); BondRunner++)
    61         if ((*BondRunner)->GetOtherAtom((*iter))->type == type2) {
     60  atom *Walker = mol->start;
     61  while (Walker->next != mol->end) {
     62    Walker = Walker->next;
     63    if (Walker->type == type1)
     64      for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++)
     65        if ((*BondRunner)->GetOtherAtom(Walker)->type == type2) {
    6266          const double distance = (*BondRunner)->GetDistanceSquared();
    6367          if (Min > distance)
     
    122126int CountHydrogenBridgeBonds(MoleculeListClass *molecules, element * InterfaceElement = NULL)
    123127{
     128  atom *Walker = NULL;
     129  atom *Runner = NULL;
    124130  int count = 0;
    125131  int OtherHydrogens = 0;
     
    127133  bool InterfaceFlag = false;
    128134  bool OtherHydrogenFlag = true;
    129   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); ++MolWalker) {
    130     molecule::iterator Walker = (*MolWalker)->begin();
    131     for(;Walker!=(*MolWalker)->end();++Walker){
    132       for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != molecules->ListOfMolecules.end(); ++MolRunner) {
    133         molecule::iterator Runner = (*MolRunner)->begin();
    134         for(;Runner!=(*MolRunner)->end();++Runner){
    135           if (((*Walker)->type->Z  == 8) && ((*Runner)->type->Z  == 8)) {
     135  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
     136    Walker = (*MolWalker)->start;
     137    while (Walker->next != (*MolWalker)->end) {
     138      Walker = Walker->next;
     139      for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != molecules->ListOfMolecules.end(); MolRunner++) {
     140        Runner = (*MolRunner)->start;
     141        while (Runner->next != (*MolRunner)->end) {
     142          Runner = Runner->next;
     143          if ((Walker->type->Z  == 8) && (Runner->type->Z  == 8)) {
    136144            // check distance
    137             const double distance = (*Runner)->x.DistanceSquared((*Walker)->x);
     145            const double distance = Runner->x.DistanceSquared(Walker->x);
    138146            if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means  different atoms
    139147              // on other atom(Runner) we check for bond to interface element and
     
    143151              OtherHydrogens = 0;
    144152              InterfaceFlag = (InterfaceElement == NULL);
    145               for (BondList::const_iterator BondRunner = (*Runner)->ListOfBonds.begin(); BondRunner != (*Runner)->ListOfBonds.end(); BondRunner++) {
    146                 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Runner);
     153              for (BondList::const_iterator BondRunner = Runner->ListOfBonds.begin(); BondRunner != Runner->ListOfBonds.end(); BondRunner++) {
     154                atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Runner);
    147155                // if hydrogen, check angle to be greater(!) than 30 degrees
    148156                if (OtherAtom->type->Z == 1) {
    149                   const double angle = CalculateAngle(&OtherAtom->x, &(*Runner)->x, &(*Walker)->x);
     157                  const double angle = CalculateAngle(&OtherAtom->x, &Runner->x, &Walker->x);
    150158                  OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON);
    151159                  Otherangle += angle;
     
    168176              if (InterfaceFlag && OtherHydrogenFlag) {
    169177                // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
    170                 for (BondList::const_iterator BondRunner = (*Walker)->ListOfBonds.begin(); BondRunner != (*Walker)->ListOfBonds.end(); BondRunner++) {
    171                   atom * const OtherAtom = (*BondRunner)->GetOtherAtom(*Walker);
     178                for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {
     179                  atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
    172180                  if (OtherAtom->type->Z == 1) {
    173181                    // check angle
    174                     if (CheckHydrogenBridgeBondAngle(*Walker, OtherAtom, *Runner)) {
    175                       DoLog(1) && (Log() << Verbose(1) << (*Walker)->getName() << ", " << OtherAtom->getName() << " and " << (*Runner)->getName() << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(&OtherAtom->x, &(*Walker)->x, &(*Runner)->x)*(180./M_PI) << "." << endl);
     182                    if (CheckHydrogenBridgeBondAngle(Walker, OtherAtom, Runner)) {
     183                      DoLog(1) && (Log() << Verbose(1) << Walker->getName() << ", " << OtherAtom->getName() << " and " << Runner->getName() << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(&OtherAtom->x, &Walker->x, &Runner->x)*(180./M_PI) << "." << endl);
    176184                      count++;
    177185                      break;
     
    197205int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second)
    198206{
     207  atom *Walker = NULL;
    199208  int count = 0;
    200209
    201210  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
    202     molecule::iterator Walker = (*MolWalker)->begin();
    203     for(;Walker!=(*MolWalker)->end();++Walker){
    204       atom * theAtom = *Walker;
    205       if ((theAtom->type == first) || (theAtom->type == second)) {  // first element matches
    206         for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
    207           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
    208           if (((OtherAtom->type == first) || (OtherAtom->type == second)) && (theAtom->nr < OtherAtom->nr)) {
     211    Walker = (*MolWalker)->start;
     212    while (Walker->next != (*MolWalker)->end) {
     213      Walker = Walker->next;
     214      if ((Walker->type == first) || (Walker->type == second)) {  // first element matches
     215        for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {
     216          atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
     217          if (((OtherAtom->type == first) || (OtherAtom->type == second)) && (Walker->nr < OtherAtom->nr)) {
    209218            count++;
    210219            DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << " bond found between " << *Walker << " and " << *OtherAtom << "." << endl);
     
    231240  bool MatchFlag[2];
    232241  bool result = false;
     242  atom *Walker = NULL;
    233243  const element * ElementArray[2];
    234244  ElementArray[0] = first;
     
    236246
    237247  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
    238     molecule::iterator Walker = (*MolWalker)->begin();
    239     for(;Walker!=(*MolWalker)->end();++Walker){
    240       atom *theAtom = *Walker;
    241       if (theAtom->type == second) {  // first element matches
     248    Walker = (*MolWalker)->start;
     249    while (Walker->next != (*MolWalker)->end) {
     250      Walker = Walker->next;
     251      if (Walker->type == second) {  // first element matches
    242252        for (int i=0;i<2;i++)
    243253          MatchFlag[i] = false;
    244         for (BondList::const_iterator BondRunner = theAtom->ListOfBonds.begin(); BondRunner != theAtom->ListOfBonds.end(); BondRunner++) {
    245           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(theAtom);
     254        for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {
     255          atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
    246256          for (int i=0;i<2;i++)
    247257            if ((!MatchFlag[i]) && (OtherAtom->type == ElementArray[i])) {
  • src/analysis_correlation.cpp

    ra7b761b r8f215d  
    4040  }
    4141  outmap = new PairCorrelationMap;
    42   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
     42  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    4343    if ((*MolWalker)->ActiveFlag) {
    4444      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    45       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    46       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    47         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    48         if ((type1 == NULL) || ((*iter)->type == type1)) {
    49           for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
     45      atom *Walker = (*MolWalker)->start;
     46      while (Walker->next != (*MolWalker)->end) {
     47        Walker = Walker->next;
     48        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     49        if ((type1 == NULL) || (Walker->type == type1)) {
     50          for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
    5051            if ((*MolOtherWalker)->ActiveFlag) {
    5152              DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
    52               for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
    53                 DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
    54                 if ((*iter)->getId() < (*runner)->getId()){
    55                   if ((type2 == NULL) || ((*runner)->type == type2)) {
    56                     distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  World::getInstance().getDomain());
    57                     //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    58                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
     53              atom *OtherWalker = (*MolOtherWalker)->start;
     54              while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
     55                OtherWalker = OtherWalker->next;
     56                DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
     57                if (Walker->nr < OtherWalker->nr)
     58                  if ((type2 == NULL) || (OtherWalker->type == type2)) {
     59                    distance = Walker->node->PeriodicDistance(*OtherWalker->node, World::getInstance().getDomain());
     60                    //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
     61                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
    5962                  }
    60                 }
    6163              }
    62             }
    6364          }
    6465        }
    6566      }
    6667    }
    67   }
     68
    6869  return outmap;
    6970};
     
    100101      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    101102      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    102       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    103         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    104         if ((type1 == NULL) || ((*iter)->type == type1)) {
    105           periodicX = *(*iter)->node;
     103      atom *Walker = (*MolWalker)->start;
     104      while (Walker->next != (*MolWalker)->end) {
     105        Walker = Walker->next;
     106        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     107        if ((type1 == NULL) || (Walker->type == type1)) {
     108          periodicX = *(Walker->node);
    106109          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    107110          // go through every range in xyz and get distance
     
    114117                  if ((*MolOtherWalker)->ActiveFlag) {
    115118                    DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
    116                     for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
    117                       DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
    118                       if ((*iter)->nr < (*runner)->nr)
    119                         if ((type2 == NULL) || ((*runner)->type == type2)) {
    120                           periodicOtherX = *(*runner)->node;
     119                    atom *OtherWalker = (*MolOtherWalker)->start;
     120                    while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
     121                      OtherWalker = OtherWalker->next;
     122                      DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl);
     123                      if (Walker->nr < OtherWalker->nr)
     124                        if ((type2 == NULL) || (OtherWalker->type == type2)) {
     125                          periodicOtherX = *(OtherWalker->node);
    121126                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    122127                          // go through every range in xyz and get distance
     
    127132                                checkOtherX.MatrixMultiplication(FullMatrix);
    128133                                distance = checkX.distance(checkOtherX);
    129                                 //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    130                                 outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
     134                                //Log() << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
     135                                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
    131136                              }
    132137                        }
     
    164169    if ((*MolWalker)->ActiveFlag) {
    165170      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    166       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    167         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    168         if ((type == NULL) || ((*iter)->type == type)) {
    169           distance = (*iter)->node->PeriodicDistance(*point, World::getInstance().getDomain());
     171      atom *Walker = (*MolWalker)->start;
     172      while (Walker->next != (*MolWalker)->end) {
     173        Walker = Walker->next;
     174        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     175        if ((type == NULL) || (Walker->type == type)) {
     176          distance = Walker->node->PeriodicDistance(*point, World::getInstance().getDomain());
    170177          DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    171           outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
     178          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
    172179        }
    173180      }
     
    204211      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    205212      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    206       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    207         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    208         if ((type == NULL) || ((*iter)->type == type)) {
    209           periodicX = *(*iter)->node;
     213      atom *Walker = (*MolWalker)->start;
     214      while (Walker->next != (*MolWalker)->end) {
     215        Walker = Walker->next;
     216        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     217        if ((type == NULL) || (Walker->type == type)) {
     218          periodicX = *(Walker->node);
    210219          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    211220          // go through every range in xyz and get distance
     
    217226                distance = checkX.distance(*point);
    218227                DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    219                 outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
     228                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
    220229              }
    221230        }
     
    252261    if ((*MolWalker)->ActiveFlag) {
    253262      DoLog(1) && (Log() << Verbose(1) << "Current molecule is " << (*MolWalker)->name << "." << endl);
    254       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    255         //Log() << Verbose(3) << "Current atom is " << *(*iter) << "." << endl;
    256         if ((type == NULL) || ((*iter)->type == type)) {
    257           TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
     263      atom *Walker = (*MolWalker)->start;
     264      while (Walker->next != (*MolWalker)->end) {
     265        Walker = Walker->next;
     266        //Log() << Verbose(1) << "Current atom is " << *Walker << "." << endl;
     267        if ((type == NULL) || (Walker->type == type)) {
     268          TriangleIntersectionList Intersections(Walker->node,Surface,LC);
    258269          distance = Intersections.GetSmallestDistance();
    259270          triangle = Intersections.GetClosestTriangle();
    260           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
     271          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
    261272        }
    262273      }
     
    303314      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    304315      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    305       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    306         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    307         if ((type == NULL) || ((*iter)->type == type)) {
    308           periodicX = *(*iter)->node;
     316      atom *Walker = (*MolWalker)->start;
     317      while (Walker->next != (*MolWalker)->end) {
     318        Walker = Walker->next;
     319        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << *Walker << "." << endl);
     320        if ((type == NULL) || (Walker->type == type)) {
     321          periodicX = *(Walker->node);
    309322          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    310323          // go through every range in xyz and get distance
     
    324337              }
    325338          // insert
    326           outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
     339          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (Walker, ShortestTriangle) ) );
    327340          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
    328341        }
  • src/atom.cpp

    ra7b761b r8f215d  
    5959atom::~atom()
    6060{
     61  unlink(this);
    6162};
    6263
  • src/bondgraph.cpp

    ra7b761b r8f215d  
    9090bool status = true;
    9191
    92   if (mol->empty()) // only construct if molecule is not empty
     92  if (mol->start->next == mol->end) // only construct if molecule is not empty
    9393    return false;
    9494
     
    124124  max_distance = 0.;
    125125
    126   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    127     if ((*iter)->type->CovalentRadius > max_distance)
    128       max_distance = (*iter)->type->CovalentRadius;
     126  atom *Runner = mol->start;
     127  while (Runner->next != mol->end) {
     128    Runner = Runner->next;
     129    if (Runner->type->CovalentRadius > max_distance)
     130      max_distance = Runner->type->CovalentRadius;
    129131  }
    130132  max_distance *= 2.;
  • src/boundary.cpp

    ra7b761b r8f215d  
    139139{
    140140        Info FunctionInfo(__func__);
     141  atom *Walker = NULL;
    141142  PointMap PointsOnBoundary;
    142143  LineMap LinesOnBoundary;
     
    164165
    165166    // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    166     for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    167       ProjectedVector = (*iter)->x - (*MolCenter);
     167    Walker = mol->start;
     168    while (Walker->next != mol->end) {
     169      Walker = Walker->next;
     170      ProjectedVector = Walker->x - (*MolCenter);
    168171      ProjectedVector.ProjectOntoPlane(AxisVector);
    169172
     
    179182        angle = 2. * M_PI - angle;
    180183      }
    181     DoLog(1) && (Log() << Verbose(1) << "Inserting " << **iter << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl);
    182       BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, (*iter))));
     184      DoLog(1) && (Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl);
     185      BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker)));
    183186      if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity
    184187        DoLog(2) && (Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl);
    185188        DoLog(2) && (Log() << Verbose(2) << "Present vector: " << *BoundaryTestPair.first->second.second << endl);
    186         DoLog(2) && (Log() << Verbose(2) << "New vector: " << **iter << endl);
     189        DoLog(2) && (Log() << Verbose(2) << "New vector: " << *Walker << endl);
    187190        const double ProjectedVectorNorm = ProjectedVector.NormSquared();
    188191        if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) {
    189192          BoundaryTestPair.first->second.first = ProjectedVectorNorm;
    190           BoundaryTestPair.first->second.second = (*iter);
     193          BoundaryTestPair.first->second.second = Walker;
    191194          DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl);
    192195        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
    193           helper = (*iter)->x;
    194           helper -= *MolCenter;
     196          helper = Walker->x - (*MolCenter);
    195197          const double oldhelperNorm = helper.NormSquared();
    196198          helper = BoundaryTestPair.first->second.second->x - (*MolCenter);
    197199          if (helper.NormSquared() < oldhelperNorm) {
    198             BoundaryTestPair.first->second.second = (*iter);
     200            BoundaryTestPair.first->second.second = Walker;
    199201            DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl);
    200202          } else {
     
    693695  int repetition[NDIM] = { 1, 1, 1 };
    694696  int TotalNoClusters = 1;
     697  atom *Walker = NULL;
    695698  double totalmass = 0.;
    696699  double clustervolume = 0.;
     
    716719
    717720  // sum up the atomic masses
    718   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    719       totalmass += (*iter)->type->mass;
     721  Walker = mol->start;
     722  while (Walker->next != mol->end) {
     723      Walker = Walker->next;
     724      totalmass += Walker->type->mass;
    720725  }
    721726  DoLog(0) && (Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl);
     
    799804  Vector Inserter;
    800805  double FillIt = false;
     806  atom *Walker = NULL;
    801807  bond *Binder = NULL;
    802808  double phi[NDIM];
     
    805811
    806812  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
    807     if ((*ListRunner)->getAtomCount() > 0) {
     813    if ((*ListRunner)->AtomCount > 0) {
    808814      DoLog(1) && (Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl);
    809815      LCList[(*ListRunner)] = new LinkedCell((*ListRunner), 10.); // get linked cell list
     
    823829  }
    824830
    825   atom * CopyAtoms[filler->getAtomCount()];
     831  filler->CountAtoms();
     832  atom * CopyAtoms[filler->AtomCount];
    826833
    827834  // calculate filler grid in [0,1]^3
     
    848855
    849856        // go through all atoms
    850         for (int i=0;i<filler->getAtomCount();i++)
     857        for (int i=0;i<filler->AtomCount;i++)
    851858          CopyAtoms[i] = NULL;
    852         for(molecule::iterator iter = filler->begin(); iter !=filler->end();++filler){
     859        Walker = filler->start;
     860        while (Walker->next != filler->end) {
     861          Walker = Walker->next;
    853862
    854863          // create atomic random translation vector ...
     
    874883
    875884          // ... and put at new position
    876           Inserter = (*iter)->x;
     885          Inserter = Walker->x;
    877886          if (DoRandomRotation)
    878887            Inserter.MatrixMultiplication(Rotations);
     
    898907            DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl);
    899908            // copy atom ...
    900             CopyAtoms[(*iter)->nr] = (*iter)->clone();
    901             CopyAtoms[(*iter)->nr]->x = Inserter;
    902             Filling->AddAtom(CopyAtoms[(*iter)->nr]);
    903             DoLog(4) && (Log() << Verbose(4) << "Filling atom " << **iter << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->nr]->x) << "." << endl);
     909            CopyAtoms[Walker->nr] = Walker->clone();
     910            CopyAtoms[Walker->nr]->x = Inserter;
     911            Filling->AddAtom(CopyAtoms[Walker->nr]);
     912            DoLog(4) && (Log() << Verbose(4) << "Filling atom " << *Walker << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[Walker->nr]->x) << "." << endl);
    904913          } else {
    905914            DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance." << endl);
    906             CopyAtoms[(*iter)->nr] = NULL;
     915            CopyAtoms[Walker->nr] = NULL;
    907916            continue;
    908917          }
     
    943952  bool TesselationFailFlag = false;
    944953
    945   mol->getAtomCount();
    946 
    947954  if (TesselStruct == NULL) {
    948955    DoLog(1) && (Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl);
     
    10161023//  // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles
    10171024//  //->InsertStraddlingPoints(mol, LCList);
    1018 //  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     1025//  mol->GoToFirst();
    10191026//  class TesselPoint *Runner = NULL;
    1020 //    Runner = *iter;
     1027//  while (!mol->IsEnd()) {
     1028//    Runner = mol->GetPoint();
    10211029//    Log() << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl;
    10221030//    if (!->IsInnerPoint(Runner, LCList)) {
     
    10261034//      Log() << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl;
    10271035//    }
     1036//    mol->GoToNext();
    10281037//  }
    10291038
     
    10341043  status = CheckListOfBaselines(TesselStruct);
    10351044
    1036   cout << "before correction" << endl;
    1037 
    10381045  // store before correction
    1039   StoreTrianglesinFile(mol, TesselStruct, filename, "");
     1046  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
    10401047
    10411048//  // correct degenerated polygons
     
    10471054  // write final envelope
    10481055  CalculateConcavityPerBoundaryPoint(TesselStruct);
    1049   cout << "after correction" << endl;
    1050   StoreTrianglesinFile(mol, TesselStruct, filename, "");
     1056  StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, "");
    10511057
    10521058  if (freeLC)
  • src/builder.cpp

    ra7b761b r8f215d  
    17441744                if (first->type != NULL) {
    17451745                  mol->AddAtom(first);  // add to molecule
    1746                   if ((configPresent == empty) && (mol->getAtomCount() != 0))
     1746                  if ((configPresent == empty) && (mol->AtomCount != 0))
    17471747                    configPresent = present;
    17481748                } else
     
    17731773                DoLog(1) && (Log() << Verbose(1) << "Depth-First-Search Analysis." << endl);
    17741774                MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
    1775                 int *MinimumRingSize = new int[mol->getAtomCount()];
     1775                int *MinimumRingSize = new int[mol->AtomCount];
    17761776                atom ***ListOfLocalAtoms = NULL;
    17771777                class StackClass<bond *> *BackEdgeStack = NULL;
     
    19211921                          int counter  = 0;
    19221922                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    1923                             if ((Boundary == NULL) || (Boundary->getAtomCount() < (*BigFinder)->getAtomCount())) {
     1923                            if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
    19241924                              Boundary = *BigFinder;
    19251925                            }
     
    19801980                performCriticalExit();
    19811981              } else {
    1982                 mol->getAtomCount();
    19831982                SaveFlag = true;
    19841983                DoLog(1) && (Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl);
     
    20992098                int counter  = 0;
    21002099                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    2101                   if ((Boundary == NULL) || (Boundary->getAtomCount() < (*BigFinder)->getAtomCount())) {
     2100                  (*BigFinder)->CountAtoms();
     2101                  if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
    21022102                    Boundary = *BigFinder;
    21032103                  }
    21042104                  counter++;
    21052105                }
    2106                 DoLog(1) && (Log() << Verbose(1) << "Biggest molecule has " << Boundary->getAtomCount() << " atoms." << endl);
     2106                DoLog(1) && (Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl);
    21072107                start = clock();
    21082108                LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.);
     
    21802180                double tmp1 = atof(argv[argptr+1]);
    21812181                atom *third = mol->FindAtom(atoi(argv[argptr]));
    2182                 if(third){
    2183                   for(molecule::iterator iter = mol->begin(); iter != mol->end();){
    2184                     if ((*iter)->x.DistanceSquared(third->x) > tmp1*tmp1){ // distance to first above radius ...
    2185                       mol->RemoveAtom(*(iter++));
    2186                     }
    2187                     else{
    2188                       ++iter;
    2189                     }
     2182                atom *first = mol->start;
     2183                if ((third != NULL) && (first != mol->end)) {
     2184                  atom *second = first->next;
     2185                  while(second != mol->end) {
     2186                    first = second;
     2187                    second = first->next;
     2188                    if (first->x.DistanceSquared(third->x) > tmp1*tmp1) // distance to first above radius ...
     2189                      mol->RemoveAtom(first);
    21902190                  }
    21912191                } else {
     
    23322332                performCriticalExit();
    23332333              } else {
    2334                 mol->getAtomCount();
    23352334                SaveFlag = true;
    23362335                DoLog(1) && (Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl);
     
    24522451                    faktor = 1;
    24532452                  }
    2454                   if (mol->getAtomCount() != 0) {  // if there is more than none
    2455                     count = mol->getAtomCount();   // is changed becausing of adding, thus has to be stored away beforehand
     2453                  mol->CountAtoms();  // recount atoms
     2454                  if (mol->AtomCount != 0) {  // if there is more than none
     2455                    count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
    24562456                    Elements = new const element *[count];
    24572457                    vectors = new Vector *[count];
    24582458                    j = 0;
    2459                     for(molecule::iterator iter = mol->begin();iter!=mol->end();++iter){
    2460                       Elements[j] = (*iter)->type;
    2461                       vectors[j] = &(*iter)->x;
     2459                    first = mol->start;
     2460                    while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
     2461                      first = first->next;
     2462                      Elements[j] = first->type;
     2463                      vectors[j] = &first->x;
    24622464                      j++;
    24632465                    }
     
    25542556int main(int argc, char **argv)
    25552557{
    2556     // while we are non interactive, we want to abort from asserts
    2557     //ASSERT_DO(Assert::Abort);
    25582558    molecule *mol = NULL;
    25592559    config *configuration = new config;
     
    25952595    }
    25962596
    2597     // In the interactive mode, we can leave the user the choice in case of error
    2598     ASSERT_DO(Assert::Ask);
    2599 
    26002597    {
    26012598      cout << ESPACKVersion << endl;
  • src/config.cpp

    ra7b761b r8f215d  
    15471547  int AtomNo = -1;
    15481548  int MolNo = 0;
     1549  atom *Walker = NULL;
    15491550  FILE *f = NULL;
    15501551
     
    15591560  fprintf(f, "# Created by MoleCuilder\n");
    15601561
    1561   for (MoleculeList::const_iterator MolRunner = MolList->ListOfMolecules.begin(); MolRunner != MolList->ListOfMolecules.end(); MolRunner++) {
     1562  for (MoleculeList::const_iterator Runner = MolList->ListOfMolecules.begin(); Runner != MolList->ListOfMolecules.end(); Runner++) {
     1563    Walker = (*Runner)->start;
    15621564    int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo");
    15631565    AtomNo = 0;
    1564     for (molecule::const_iterator iter = (*MolRunner)->begin(); iter != (*MolRunner)->end(); ++iter) {
    1565       sprintf(name, "%2s%2d",(*iter)->type->symbol, elementNo[(*iter)->type->Z]);
    1566       elementNo[(*iter)->type->Z] = (elementNo[(*iter)->type->Z]+1) % 100;   // confine to two digits
     1566    while (Walker->next != (*Runner)->end) {
     1567      Walker = Walker->next;
     1568      sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
     1569      elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
    15671570      fprintf(f,
    15681571             "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
    1569              (*iter)->nr,                /* atom serial number */
     1572             Walker->nr,                /* atom serial number */
    15701573             name,         /* atom name */
    1571              (*MolRunner)->name,      /* residue name */
     1574             (*Runner)->name,      /* residue name */
    15721575             'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
    15731576             MolNo,         /* residue sequence number */
    1574              (*iter)->node->at(0),                 /* position X in Angstroem */
    1575              (*iter)->node->at(1),                 /* position Y in Angstroem */
    1576              (*iter)->node->at(2),                 /* position Z in Angstroem */
    1577              (double)(*iter)->type->Valence,         /* occupancy */
    1578              (double)(*iter)->type->NoValenceOrbitals,          /* temperature factor */
     1577             Walker->node->at(0),                 /* position X in Angstroem */
     1578             Walker->node->at(1),                 /* position Y in Angstroem */
     1579             Walker->node->at(2),                 /* position Z in Angstroem */
     1580             (double)Walker->type->Valence,         /* occupancy */
     1581             (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
    15791582             "0",            /* segment identifier */
    1580              (*iter)->type->symbol,    /* element symbol */
     1583             Walker->type->symbol,    /* element symbol */
    15811584             "0");           /* charge */
    15821585      AtomNo++;
     
    15981601{
    15991602  int AtomNo = -1;
     1603  atom *Walker = NULL;
    16001604  FILE *f = NULL;
    16011605
     
    16121616  fprintf(f, "# Created by MoleCuilder\n");
    16131617
     1618  Walker = mol->start;
    16141619  AtomNo = 0;
    1615   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    1616     sprintf(name, "%2s%2d",(*iter)->type->symbol, elementNo[(*iter)->type->Z]);
    1617     elementNo[(*iter)->type->Z] = (elementNo[(*iter)->type->Z]+1) % 100;   // confine to two digits
     1620  while (Walker->next != mol->end) {
     1621    Walker = Walker->next;
     1622    sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
     1623    elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
    16181624    fprintf(f,
    16191625           "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
    1620            (*iter)->nr,                /* atom serial number */
     1626           Walker->nr,                /* atom serial number */
    16211627           name,         /* atom name */
    16221628           mol->name,      /* residue name */
    16231629           'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
    16241630           0,         /* residue sequence number */
    1625            (*iter)->node->at(0),                 /* position X in Angstroem */
    1626            (*iter)->node->at(1),                 /* position Y in Angstroem */
    1627            (*iter)->node->at(2),                 /* position Z in Angstroem */
    1628            (double)(*iter)->type->Valence,         /* occupancy */
    1629            (double)(*iter)->type->NoValenceOrbitals,          /* temperature factor */
     1631           Walker->node->at(0),                 /* position X in Angstroem */
     1632           Walker->node->at(1),                 /* position Y in Angstroem */
     1633           Walker->node->at(2),                 /* position Z in Angstroem */
     1634           (double)Walker->type->Valence,         /* occupancy */
     1635           (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
    16301636           "0",            /* segment identifier */
    1631            (*iter)->type->symbol,    /* element symbol */
     1637           Walker->type->symbol,    /* element symbol */
    16321638           "0");           /* charge */
    16331639    AtomNo++;
     
    16471653bool config::SaveTREMOLO(const char * const filename, const molecule * const mol) const
    16481654{
     1655  atom *Walker = NULL;
    16491656  ofstream *output = NULL;
    16501657  stringstream * const fname = new stringstream;
     
    16591666
    16601667  // scan maximum number of neighbours
     1668  Walker = mol->start;
    16611669  int MaxNeighbours = 0;
    1662   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    1663     const int count = (*iter)->ListOfBonds.size();
     1670  while (Walker->next != mol->end) {
     1671    Walker = Walker->next;
     1672    const int count = Walker->ListOfBonds.size();
    16641673    if (MaxNeighbours < count)
    16651674      MaxNeighbours = count;
    16661675  }
    1667   *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;
    1668 
    1669   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    1670     *output << (*iter)->nr << "\t";
    1671     *output << (*iter)->getName() << "\t";
     1676  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
     1677
     1678  Walker = mol->start;
     1679  while (Walker->next != mol->end) {
     1680    Walker = Walker->next;
     1681    *output << Walker->nr << "\t";
     1682    *output << Walker->getName() << "\t";
    16721683    *output << mol->name << "\t";
    16731684    *output << 0 << "\t";
    1674     *output << (*iter)->node->at(0) << "\t" << (*iter)->node->at(1) << "\t" << (*iter)->node->at(2) << "\t";
    1675     *output << static_cast<double>((*iter)->type->Valence) << "\t";
    1676     *output << (*iter)->type->symbol << "\t";
    1677     for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
    1678       *output << (*runner)->GetOtherAtom(*iter)->nr << "\t";
    1679     for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
     1685    *output << Walker->node->at(0) << "\t" << Walker->node->at(1) << "\t" << Walker->node->at(2) << "\t";
     1686    *output << static_cast<double>(Walker->type->Valence) << "\t";
     1687    *output << Walker->type->symbol << "\t";
     1688    for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
     1689      *output << (*runner)->GetOtherAtom(Walker)->nr << "\t";
     1690    for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
    16801691      *output << "-\t";
    16811692    *output << endl;
     
    16971708bool config::SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const
    16981709{
     1710  atom *Walker = NULL;
    16991711  ofstream *output = NULL;
    17001712  stringstream * const fname = new stringstream;
     
    17111723  int MaxNeighbours = 0;
    17121724  for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
    1713     for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    1714       const int count = (*iter)->ListOfBonds.size();
     1725    Walker = (*MolWalker)->start;
     1726    while (Walker->next != (*MolWalker)->end) {
     1727      Walker = Walker->next;
     1728      const int count = Walker->ListOfBonds.size();
    17151729      if (MaxNeighbours < count)
    17161730        MaxNeighbours = count;
    17171731    }
    17181732  }
    1719   *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;
     1733  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
    17201734
    17211735  // create global to local id map
     
    17381752    int AtomNo = 0;
    17391753    for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
    1740       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
     1754      Walker = (*MolWalker)->start;
     1755      while (Walker->next != (*MolWalker)->end) {
     1756        Walker = Walker->next;
    17411757        *output << AtomNo+1 << "\t";
    1742         *output << (*iter)->getName() << "\t";
     1758        *output << Walker->getName() << "\t";
    17431759        *output << (*MolWalker)->name << "\t";
    17441760        *output << MolCounter+1 << "\t";
    1745         *output << (*iter)->node->at(0) << "\t" << (*iter)->node->at(1) << "\t" << (*iter)->node->at(2) << "\t";
    1746         *output << (double)(*iter)->type->Valence << "\t";
    1747         *output << (*iter)->type->symbol << "\t";
    1748         for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
    1749           *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom(*iter)->nr ]+1 << "\t";
    1750         for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
     1761        *output << Walker->node->at(0) << "\t" << Walker->node->at(1) << "\t" << Walker->node->at(2) << "\t";
     1762        *output << (double)Walker->type->Valence << "\t";
     1763        *output << Walker->type->symbol << "\t";
     1764        for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
     1765          *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom(Walker)->nr ]+1 << "\t";
     1766        for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
    17511767          *output << "-\t";
    17521768        *output << endl;
  • src/helpers.hpp

    ra7b761b r8f215d  
    8383};
    8484
     85/** Creates a lookup table for true father's Atom::Nr -> atom ptr.
     86 * \param *start begin of chain list
     87 * \paran *end end of chain list
     88 * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
     89 * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
     90 * \return true - success, false - failure
     91 */
     92template <typename T> bool CreateFatherLookupTable(T *start, T *end, T **&LookupTable, int count = 0)
     93{
     94  bool status = true;
     95  T *Walker = NULL;
     96  int AtomNo;
     97
     98  if (LookupTable != NULL) {
     99    DoLog(0) && (Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl);
     100    return false;
     101  }
     102
     103  // count them
     104  if (count == 0) {
     105    Walker = start;
     106    while (Walker->next != end) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
     107      Walker = Walker->next;
     108      count = (count < Walker->GetTrueFather()->nr) ? Walker->GetTrueFather()->nr : count;
     109    }
     110  }
     111  if (count <= 0) {
     112    DoLog(0) && (Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl);
     113    return false;
     114  }
     115
     116  // allocate and fill
     117  LookupTable = Calloc<T*>(count, "CreateFatherLookupTable - **LookupTable");
     118  if (LookupTable == NULL) {
     119    DoeLog(0) && (eLog()<< Verbose(0) << "LookupTable memory allocation failed!" << endl);
     120    performCriticalExit();
     121    status = false;
     122  } else {
     123    Walker = start;
     124    while (Walker->next != end) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
     125      Walker = Walker->next;
     126      AtomNo = Walker->GetTrueFather()->nr;
     127      if ((AtomNo >= 0) && (AtomNo < count)) {
     128        //*out << "Setting LookupTable[" << AtomNo << "] to " << *Walker << endl;
     129        LookupTable[AtomNo] = Walker;
     130      } else {
     131        DoLog(0) && (Log() << Verbose(0) << "Walker " << *Walker << " exceeded range of nuclear ids [0, " << count << ")." << endl);
     132        status = false;
     133        break;
     134      }
     135    }
     136  }
     137
     138  return status;
     139};
     140
    85141/** Frees a two-dimensional array.
    86142 * \param *ptr pointer to array
  • src/lists.hpp

    ra7b761b r8f215d  
    134134};
    135135
     136/** Returns the first marker in a chain list.
     137 * \param *me one arbitrary item in chain list
     138 * \return poiner to first marker
     139 */
     140template <typename X> X *GetFirst(X *me)
     141{
     142  X *Binder = me;
     143  while(Binder->previous != 0)
     144    Binder = Binder->previous;
     145  return Binder;
     146};
     147
     148/** Returns the last marker in a chain list.
     149 * \param *me one arbitrary item in chain list
     150 * \return poiner to last marker
     151 */
     152template <typename X> X *GetLast(X *me)
     153{
     154  X *Binder = me;
     155  while(Binder->next != 0)
     156    Binder = Binder->next;
     157  return Binder;
     158};
     159
    136160#endif /* LISTS_HPP_ */
  • src/molecule.cpp

    ra7b761b r8f215d  
    3535 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
    3636 */
    37 molecule::molecule(const periodentafel * const teil) : elemente(teil),
    38   first(new bond(0, 0, 1, -1)), last(new bond(0, 0, 1, -1)), MDSteps(0),
     37molecule::molecule(const periodentafel * const teil) : elemente(teil), start(World::getInstance().createAtom()), end(World::getInstance().createAtom()),
     38  first(new bond(start, end, 1, -1)), last(new bond(start, end, 1, -1)), MDSteps(0), AtomCount(0),
    3939  BondCount(0), ElementCount(0), NoNonHydrogen(0), NoNonBonds(0), NoCyclicBonds(0), BondDistance(0.),
    4040  ActiveFlag(false), IndexNr(-1),
    4141  formula(this,boost::bind(&molecule::calcFormula,this)),
    42   AtomCount(this,boost::bind(&molecule::doCountAtoms,this)), last_atom(0),  InternalPointer(begin())
    43 {
     42  last_atom(0),
     43  InternalPointer(start)
     44{
     45  // init atom chain list
     46  start->father = NULL;
     47  end->father = NULL;
     48  link(start,end);
     49
    4450  // init bond chain list
    4551  link(first,last);
     
    6369  delete(first);
    6470  delete(last);
     71  end->getWorld()->destroyAtom(end);
     72  start->getWorld()->destroyAtom(start);
    6573};
    6674
     
    7583}
    7684
    77 int molecule::getAtomCount() const{
    78   return *AtomCount;
    79 }
    80 
    8185void molecule::setName(const std::string _name){
    8286  OBSERVE;
     
    100104  stringstream sstr;
    101105  periodentafel *periode = World::getInstance().getPeriode();
    102   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    103     counts[(*iter)->type->getNumber()]++;
     106  for(atom *Walker = start; Walker != end; Walker = Walker->next) {
     107    counts[Walker->type->getNumber()]++;
    104108  }
    105109  std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
     
    111115}
    112116
    113 /************************** Access to the List of Atoms ****************/
    114 
    115 
    116 molecule::iterator molecule::begin(){
    117   return molecule::iterator(atoms.begin(),this);
    118 }
    119 
    120 molecule::const_iterator molecule::begin() const{
    121   return atoms.begin();
    122 }
    123 
    124 molecule::iterator molecule::end(){
    125   return molecule::iterator(atoms.end(),this);
    126 }
    127 
    128 molecule::const_iterator molecule::end() const{
    129   return atoms.end();
    130 }
    131 
    132 bool molecule::empty() const
    133 {
    134   return (begin() == end());
    135 }
    136 
    137 size_t molecule::size() const
    138 {
    139   size_t counter = 0;
    140   for (molecule::const_iterator iter = begin(); iter != end (); ++iter)
    141     counter++;
    142   return counter;
    143 }
    144 
    145 molecule::const_iterator molecule::erase( const_iterator loc )
    146 {
    147   molecule::const_iterator iter = loc;
    148   iter--;
    149   atoms.erase( loc );
    150   return iter;
    151 }
    152 
    153 molecule::const_iterator molecule::erase( atom *& key )
    154 {
    155   cout << "trying to erase atom" << endl;
    156   molecule::const_iterator iter = find(key);
    157   if (iter != end()){
    158     // remove this position and step forward (post-increment)
    159     atoms.erase( iter++ );
    160   }
    161   return iter;
    162 }
    163 
    164 molecule::const_iterator molecule::find ( atom *& key ) const
    165 {
    166   return atoms.find( key );
    167 }
    168 
    169 pair<molecule::iterator,bool> molecule::insert ( atom * const key )
    170 {
    171   pair<atomSet::iterator,bool> res = atoms.insert(key);
    172   return pair<iterator,bool>(iterator(res.first,this),res.second);
    173 }
    174117
    175118/** Adds given atom \a *pointer from molecule list.
     
    180123bool molecule::AddAtom(atom *pointer)
    181124{
     125  bool retval = false;
    182126  OBSERVE;
    183127  if (pointer != NULL) {
    184128    pointer->sort = &pointer->nr;
     129    pointer->nr = last_atom++;  // increase number within molecule
     130    AtomCount++;
    185131    if (pointer->type != NULL) {
    186132      if (ElementsInMolecule[pointer->type->Z] == 0)
     
    195141      }
    196142    }
    197     insert(pointer);
    198   }
    199   return true;
     143    retval = add(pointer, end);
     144  }
     145  return retval;
    200146};
    201147
     
    211157  if (pointer != NULL) {
    212158    atom *walker = pointer->clone();
    213     walker->setName(pointer->getName());
     159    stringstream sstr;
     160    sstr << pointer->getName();
     161    walker->setName(sstr.str());
    214162    walker->nr = last_atom++;  // increase number within molecule
    215     insert(walker);
     163    add(walker, end);
    216164    if ((pointer->type != NULL) && (pointer->type->Z != 1))
    217165      NoNonHydrogen++;
     166    AtomCount++;
    218167    retval=walker;
    219168  }
     
    625574
    626575  // copy values
     576  copy->CountAtoms();
    627577  copy->CountElements();
    628578  if (first->next != last) {  // if adjaceny list is present
     
    659609{
    660610  bond *Binder = NULL;
    661 
    662   // some checks to make sure we are able to create the bond
    663   ASSERT(atom1, "First atom in bond-creation was an invalid pointer");
    664   ASSERT(atom2, "Second atom in bond-creation was an invalid pointer");
    665   ASSERT(FindAtom(atom1->nr),"First atom in bond-creation was not part of molecule");
    666   ASSERT(FindAtom(atom2->nr),"Second atom in bond-creation was not parto of molecule");
    667 
    668   Binder = new bond(atom1, atom2, degree, BondCount++);
    669   atom1->RegisterBond(Binder);
    670   atom2->RegisterBond(Binder);
    671   if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
    672     NoNonBonds++;
    673   add(Binder, last);
    674 
     611  if ((atom1 != NULL) && (FindAtom(atom1->nr) != NULL) && (atom2 != NULL) && (FindAtom(atom2->nr) != NULL)) {
     612    Binder = new bond(atom1, atom2, degree, BondCount++);
     613    atom1->RegisterBond(Binder);
     614    atom2->RegisterBond(Binder);
     615    if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
     616      NoNonBonds++;
     617    add(Binder, last);
     618  } else {
     619    DoeLog(1) && (eLog()<< Verbose(1) << "Could not add bond between " << atom1->getName() << " and " << atom2->getName() << " as one or both are not present in the molecule." << endl);
     620  }
    675621  return Binder;
    676622};
     
    746692bool molecule::RemoveAtom(atom *pointer)
    747693{
    748   ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");
    749   OBSERVE;
    750694  if (ElementsInMolecule[pointer->type->Z] != 0)  { // this would indicate an error
    751695    ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
     696    AtomCount--;
    752697  } else
    753698    DoeLog(1) && (eLog()<< Verbose(1) << "Atom " << pointer->getName() << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl);
     
    755700    ElementCount--;
    756701  RemoveBonds(pointer);
    757   erase(pointer);
    758   return true;
     702  return remove(pointer, start, end);
    759703};
    760704
     
    773717  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    774718    ElementCount--;
    775   erase(pointer);
     719  unlink(pointer);
    776720  return true;
    777721};
     
    782726bool molecule::CleanupMolecule()
    783727{
    784   for (molecule::iterator iter = begin(); !empty(); iter = begin())
    785       erase(iter);
    786   return (cleanup(first,last));
     728  return (cleanup(first,last) && cleanup(start,end));
    787729};
    788730
     
    791733 * \return pointer to atom or NULL
    792734 */
    793 atom * molecule::FindAtom(int Nr)  const
    794 {
    795   molecule::const_iterator iter = begin();
    796   for (; iter != end(); ++iter)
    797     if ((*iter)->nr == Nr)
    798       break;
    799   if (iter != end()) {
     735atom * molecule::FindAtom(int Nr)  const{
     736  atom * walker = find(&Nr, start,end);
     737  if (walker != NULL) {
    800738    //Log() << Verbose(0) << "Found Atom Nr. " << walker->nr << endl;
    801     return (*iter);
     739    return walker;
    802740  } else {
    803741    DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);
     
    929867    now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    930868    for (int step=0;step<MDSteps;step++) {
    931       *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
     869      *output << AtomCount << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
    932870      ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
    933871    }
     
    946884  if (output != NULL) {
    947885    now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
    948     *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
     886    *output << AtomCount << "\n\tCreated by molecuilder on " << ctime(&now);
    949887    ActOnAllAtoms( &atom::OutputXYZLine, output );
    950888    return true;
     
    956894 * \param *out output stream for debugging
    957895 */
    958 int molecule::doCountAtoms()
    959 {
    960   int res = size();
     896void molecule::CountAtoms()
     897{
    961898  int i = 0;
    962   NoNonHydrogen = 0;
    963   for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    964     (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
    965     if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
    966       NoNonHydrogen++;
    967     stringstream sstr;
    968     sstr << (*iter)->type->symbol << (*iter)->nr+1;
    969     (*iter)->setName(sstr.str());
    970     Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl;
     899  atom *Walker = start;
     900  while (Walker->next != end) {
     901    Walker = Walker->next;
    971902    i++;
    972903  }
    973   return res;
     904  if ((AtomCount == 0) || (i != AtomCount)) {
     905    DoLog(3) && (Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl);
     906    AtomCount = i;
     907
     908    // count NonHydrogen atoms and give each atom a unique name
     909    if (AtomCount != 0) {
     910      i=0;
     911      NoNonHydrogen = 0;
     912      Walker = start;
     913      while (Walker->next != end) {
     914        Walker = Walker->next;
     915        Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
     916        if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
     917          NoNonHydrogen++;
     918        stringstream sstr;
     919        sstr << Walker->type->symbol << Walker->nr+1;
     920        Walker->setName(sstr.str());
     921        DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->getName() << "." << endl);
     922        i++;
     923      }
     924    } else
     925      DoLog(3) && (Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl);
     926  }
    974927};
    975928
     
    1033986  /// first count both their atoms and elements and update lists thereby ...
    1034987  //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
     988  CountAtoms();
     989  OtherMolecule->CountAtoms();
    1035990  CountElements();
    1036991  OtherMolecule->CountElements();
     
    1039994  /// -# AtomCount
    1040995  if (result) {
    1041     if (getAtomCount() != OtherMolecule->getAtomCount()) {
    1042       DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl);
     996    if (AtomCount != OtherMolecule->AtomCount) {
     997      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl);
    1043998      result = false;
    1044     } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
     999    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
    10451000  }
    10461001  /// -# ElementCount
     
    10791034  if (result) {
    10801035    DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
    1081     Distances = Calloc<double>(getAtomCount(), "molecule::IsEqualToWithinThreshold: Distances");
    1082     OtherDistances = Calloc<double>(getAtomCount(), "molecule::IsEqualToWithinThreshold: OtherDistances");
     1036    Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
     1037    OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
    10831038    SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
    10841039    SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
     
    10861041    /// ... sort each list (using heapsort (o(N log N)) from GSL)
    10871042    DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
    1088     PermMap = Calloc<size_t>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *PermMap");
    1089     OtherPermMap = Calloc<size_t>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *OtherPermMap");
    1090     gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);
    1091     gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);
    1092     PermutationMap = Calloc<int>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *PermutationMap");
     1043    PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
     1044    OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
     1045    gsl_heapsort_index (PermMap, Distances, AtomCount, sizeof(double), CompareDoubles);
     1046    gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
     1047    PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
    10931048    DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
    1094     for(int i=getAtomCount();i--;)
     1049    for(int i=AtomCount;i--;)
    10951050      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
    10961051
     
    10981053    DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
    10991054    flag = 0;
    1100     for (int i=0;i<getAtomCount();i++) {
     1055    for (int i=0;i<AtomCount;i++) {
    11011056      DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl);
    11021057      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
     
    11341089int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule)
    11351090{
     1091  atom *Walker = NULL, *OtherWalker = NULL;
    11361092  DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);
    1137   int *AtomicMap = Malloc<int>(getAtomCount(), "molecule::GetAtomicMap: *AtomicMap");
    1138   for (int i=getAtomCount();i--;)
     1093  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
     1094  for (int i=AtomCount;i--;)
    11391095    AtomicMap[i] = -1;
    11401096  if (OtherMolecule == this) {  // same molecule
    1141     for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
     1097    for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
    11421098      AtomicMap[i] = i;
    11431099    DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl);
    11441100  } else {
    11451101    DoLog(4) && (Log() << Verbose(4) << "Map is ");
    1146     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    1147       if ((*iter)->father == NULL) {
    1148         AtomicMap[(*iter)->nr] = -2;
     1102    Walker = start;
     1103    while (Walker->next != end) {
     1104      Walker = Walker->next;
     1105      if (Walker->father == NULL) {
     1106        AtomicMap[Walker->nr] = -2;
    11491107      } else {
    1150         for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) {
     1108        OtherWalker = OtherMolecule->start;
     1109        while (OtherWalker->next != OtherMolecule->end) {
     1110          OtherWalker = OtherWalker->next;
    11511111      //for (int i=0;i<AtomCount;i++) { // search atom
    11521112        //for (int j=0;j<OtherMolecule->AtomCount;j++) {
    1153           //Log() << Verbose(4) << "Comparing father " << (*iter)->father << " with the other one " << (*runner)->father << "." << endl;
    1154           if ((*iter)->father == (*runner))
    1155             AtomicMap[(*iter)->nr] = (*runner)->nr;
     1113          //Log() << Verbose(4) << "Comparing father " << Walker->father << " with the other one " << OtherWalker->father << "." << endl;
     1114          if (Walker->father == OtherWalker)
     1115            AtomicMap[Walker->nr] = OtherWalker->nr;
    11561116        }
    11571117      }
    1158       DoLog(0) && (Log() << Verbose(0) << AtomicMap[(*iter)->nr] << "\t");
     1118      DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t");
    11591119    }
    11601120    DoLog(0) && (Log() << Verbose(0) << endl);
     
    11901150void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const
    11911151{
    1192   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    1193     array[((*iter)->*index)] = (*iter);
     1152  atom *Walker = start;
     1153  while (Walker->next != end) {
     1154    Walker = Walker->next;
     1155    array[(Walker->*index)] = Walker;
    11941156  }
    11951157};
  • src/molecule.hpp

    ra7b761b r8f215d  
    3434#include "tesselation.hpp"
    3535#include "Patterns/Observer.hpp"
    36 #include "Patterns/ObservedIterator.hpp"
    3736#include "Patterns/Cacheable.hpp"
    3837
     
    8988  friend molecule *NewMolecule();
    9089  friend void DeleteMolecule(molecule *);
    91 
    9290  public:
    93     typedef std::set<atom*> atomSet;
    94     typedef ObservedIterator<atomSet> iterator;
    95     typedef atomSet::const_iterator const_iterator;
    96 
    9791    const periodentafel * const elemente; //!< periodic table with each element
    98     // old deprecated atom handling
    99     //atom *start;        //!< start of atom list
    100     //atom *end;          //!< end of atom list
     92    atom *start;        //!< start of atom list
     93    atom *end;          //!< end of atom list
    10194    bond *first;        //!< start of bond list
    10295    bond *last;         //!< end of bond list
    10396    int MDSteps;        //!< The number of MD steps in Trajectories
    104     //int AtomCount;          //!< number of atoms, brought up-to-date by CountAtoms()
     97    int AtomCount;          //!< number of atoms, brought up-to-date by CountAtoms()
    10598    int BondCount;          //!< number of atoms, brought up-to-date by CountBonds()
    10699    int ElementCount;       //!< how many unique elements are therein
     
    117110  private:
    118111    Cacheable<string> formula;
    119     Cacheable<int>    AtomCount;
    120112    moleculeId_t id;
    121     atomSet atoms; //<!set of atoms
    122113  protected:
    123     //void CountAtoms();
    124     /**
    125      * this iterator type should be used for internal variables, \
    126      * since it will not lock
    127      */
    128     typedef atomSet::iterator internal_iterator;
    129 
    130 
    131114    molecule(const periodentafel * const teil);
    132115    virtual ~molecule();
     
    136119  //getter and setter
    137120  const std::string getName();
    138   int getAtomCount() const;
    139   int doCountAtoms();
    140121  moleculeId_t getId();
    141122  void setId(moleculeId_t);
     
    144125  std::string calcFormula();
    145126
    146   iterator begin();
    147   const_iterator begin() const;
    148   iterator end();
    149   const_iterator end() const;
    150   bool empty() const;
    151   size_t size() const;
    152   const_iterator erase( const_iterator loc );
    153   const_iterator erase( atom *& key );
    154   const_iterator find (  atom *& key ) const;
    155   pair<iterator,bool> insert ( atom * const key );
    156 
    157127
    158128  // re-definition of virtual functions from PointCloud
     
    160130  Vector *GetCenter() const ;
    161131  TesselPoint *GetPoint() const ;
     132  TesselPoint *GetTerminalPoint() const ;
    162133  int GetMaxId() const;
    163134  void GoToNext() const ;
     135  void GoToPrevious() const ;
    164136  void GoToFirst() const ;
     137  void GoToLast() const ;
    165138  bool IsEmpty() const ;
    166139  bool IsEnd() const ;
     
    257230
    258231  /// Count and change present atoms' coordination.
     232  void CountAtoms();
    259233  void CountElements();
    260234  void CalculateOrbitals(class config &configuration);
     
    325299  bool StoreForcesFile(MoleculeListClass *BondFragments, char *path, int *SortIndex);
    326300  bool CreateMappingLabelsToConfigSequence(int *&SortIndex);
    327   bool CreateFatherLookupTable(atom **&LookupTable, int count = 0);
    328301  void BreadthFirstSearchAdd(molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem);
    329302  /// -# BOSSANOVA
     
    354327  private:
    355328  int last_atom;      //!< number given to last atom
    356   mutable internal_iterator InternalPointer;  //!< internal pointer for PointCloud
     329  mutable atom *InternalPointer;  //!< internal pointer for PointCloud
    357330};
    358331
  • src/molecule_dynamics.cpp

    ra7b761b r8f215d  
    2828  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
    2929  gsl_vector *x = gsl_vector_alloc(NDIM);
     30  atom * Runner = mol->start;
    3031  atom *Sprinter = NULL;
    3132  Vector trajectory1, trajectory2, normal, TestVector;
    3233  double Norm1, Norm2, tmp, result = 0.;
    3334
    34   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    35     if ((*iter) == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
     35  while (Runner->next != mol->end) {
     36    Runner = Runner->next;
     37    if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
    3638      break;
    3739    // determine normalized trajectories direction vector (n1, n2)
     
    4042    trajectory1.Normalize();
    4143    Norm1 = trajectory1.Norm();
    42     Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
    43     trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
     44    Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
     45    trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
    4446    trajectory2.Normalize();
    4547    Norm2 = trajectory1.Norm();
    4648    // check whether either is zero()
    4749    if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    48       tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
     50      tmp = Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.startstep));
    4951    } else if (Norm1 < MYEPSILON) {
    5052      Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
    51       trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
     53      trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
    5254      trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    5355      trajectory1 -= trajectory2;   // project the part in norm direction away
    5456      tmp = trajectory1.Norm();  // remaining norm is distance
    5557    } else if (Norm2 < MYEPSILON) {
    56       Sprinter = Params.PermutationMap[(*iter)->nr];   // find second target point
     58      Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
    5759      trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);  // copy second offset
    5860      trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
     
    6466  //        Log() << Verbose(0) << " and ";
    6567  //        Log() << Verbose(0) << trajectory2;
    66       tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
     68      tmp = Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.startstep));
    6769  //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    6870    } else { // determine distance by finding minimum distance
    69   //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear independent ";
     71  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
    7072  //        Log() << Verbose(0) << endl;
    7173  //        Log() << Verbose(0) << "First Trajectory: ";
     
    8385        gsl_matrix_set(A, 1, i, trajectory2[i]);
    8486        gsl_matrix_set(A, 2, i, normal[i]);
    85         gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - (*iter)->Trajectory.R.at(Params.startstep)[i]));
     87        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - Runner->Trajectory.R.at(Params.startstep)[i]));
    8688      }
    8789      // solve the linear system by Householder transformations
     
    9496      trajectory2.Scale(gsl_vector_get(x,1));
    9597      normal.Scale(gsl_vector_get(x,2));
    96       TestVector = (*iter)->Trajectory.R.at(Params.startstep) + trajectory2 + normal
     98      TestVector = Runner->Trajectory.R.at(Params.startstep) + trajectory2 + normal
    9799                   - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
    98100      if (TestVector.Norm() < MYEPSILON) {
     
    123125{
    124126  double result = 0.;
    125   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    126     if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[(*iter)->nr]) && (Walker->nr < (*iter)->nr)) {
     127  atom * Runner = mol->start;
     128  while (Runner->next != mol->end) {
     129    Runner = Runner->next;
     130    if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
    127131  //    atom *Sprinter = PermutationMap[Walker->nr];
    128   //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
     132  //        Log() << Verbose(0) << *Walker << " and " << *Runner << " are heading to the same target at ";
    129133  //        Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep);
    130134  //        Log() << Verbose(0) << ", penalting." << endl;
     
    157161  // go through every atom
    158162  atom *Runner = NULL;
    159   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     163  atom *Walker = start;
     164  while (Walker->next != end) {
     165    Walker = Walker->next;
    160166    // first term: distance to target
    161     Runner = Params.PermutationMap[(*iter)->nr];   // find target point
    162     tmp = ((*iter)->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
     167    Runner = Params.PermutationMap[Walker->nr];   // find target point
     168    tmp = (Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
    163169    tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    164170    result += Params.PenaltyConstants[0] * tmp;
     
    166172
    167173    // second term: sum of distances to other trajectories
    168     result += SumDistanceOfTrajectories((*iter), this, Params);
     174    result += SumDistanceOfTrajectories(Walker, this, Params);
    169175
    170176    // third term: penalty for equal targets
    171     result += PenalizeEqualTargets((*iter), this, Params);
     177    result += PenalizeEqualTargets(Walker, this, Params);
    172178  }
    173179
     
    207213void FillDistanceList(molecule *mol, struct EvaluatePotential &Params)
    208214{
    209   for (int i=mol->getAtomCount(); i--;) {
     215  for (int i=mol->AtomCount; i--;) {
    210216    Params.DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
    211217    Params.DistanceList[i]->clear();
    212218  }
    213219
    214   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    215     for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
    216       Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->Trajectory.R.at(Params.startstep).distance((*runner)->Trajectory.R.at(Params.endstep)), (*runner)) );
     220  atom *Runner = NULL;
     221  atom *Walker = mol->start;
     222  while (Walker->next != mol->end) {
     223    Walker = Walker->next;
     224    Runner = mol->start;
     225    while(Runner->next != mol->end) {
     226      Runner = Runner->next;
     227      Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)), Runner) );
    217228    }
    218229  }
     
    226237void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params)
    227238{
    228   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    229     Params.StepList[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin();    // stores the step to the next iterator that could be a possible next target
    230     Params.PermutationMap[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin()->second;   // always pick target with the smallest distance
    231     Params.DoubleList[Params.DistanceList[(*iter)->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
    232     Params.DistanceIterators[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin();    // and remember which one we picked
    233     DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << Params.DistanceList[(*iter)->nr]->begin()->first << "." << endl);
     239  atom *Walker = mol->start;
     240  while (Walker->next != mol->end) {
     241    Walker = Walker->next;
     242    Params.StepList[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
     243    Params.PermutationMap[Walker->nr] = Params.DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
     244    Params.DoubleList[Params.DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
     245    Params.DistanceIterators[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // and remember which one we picked
     246    DoLog(2) && (Log() << Verbose(2) << *Walker << " starts with distance " << Params.DistanceList[Walker->nr]->begin()->first << "." << endl);
    234247  }
    235248};
     
    272285void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params)
    273286{
    274   molecule::const_iterator iter = mol->begin();
     287  atom *Walker = mol->start;
    275288  DistanceMap::iterator NewBase;
    276289  double Potential = fabs(mol->ConstrainedPotential(Params));
    277290
    278   if (mol->empty()) {
    279     eLog() << Verbose(1) << "Molecule is empty." << endl;
    280     return;
    281   }
    282291  while ((Potential) > Params.PenaltyConstants[2]) {
    283     PrintPermutationMap(mol->getAtomCount(), Params);
    284     iter++;
    285     if (iter == mol->end()) // round-robin at the end
    286       iter = mol->begin();
    287     if (Params.DoubleList[Params.DistanceIterators[(*iter)->nr]->second->nr] <= 1)  // no need to make those injective that aren't
     292    PrintPermutationMap(mol->AtomCount, Params);
     293    Walker = Walker->next;
     294    if (Walker == mol->end) // round-robin at the end
     295      Walker = mol->start->next;
     296    if (Params.DoubleList[Params.DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
    288297      continue;
    289298    // now, try finding a new one
    290     Potential = TryNextNearestNeighbourForInjectivePermutation(mol, (*iter), Potential, Params);
    291   }
    292   for (int i=mol->getAtomCount(); i--;) // now each single entry in the DoubleList should be <=1
     299    Potential = TryNextNearestNeighbourForInjectivePermutation(mol, Walker, Potential, Params);
     300  }
     301  for (int i=mol->AtomCount; i--;) // now each single entry in the DoubleList should be <=1
    293302    if (Params.DoubleList[i] > 1) {
    294303      DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
     
    329338  double Potential, OldPotential, OlderPotential;
    330339  struct EvaluatePotential Params;
    331   Params.PermutationMap = Calloc<atom*>(getAtomCount(), "molecule::MinimiseConstrainedPotential: Params.**PermutationMap");
    332   Params.DistanceList = Malloc<DistanceMap*>(getAtomCount(), "molecule::MinimiseConstrainedPotential: Params.**DistanceList");
    333   Params.DistanceIterators = Malloc<DistanceMap::iterator>(getAtomCount(), "molecule::MinimiseConstrainedPotential: Params.*DistanceIterators");
    334   Params.DoubleList = Calloc<int>(getAtomCount(), "molecule::MinimiseConstrainedPotential: Params.*DoubleList");
    335   Params.StepList = Malloc<DistanceMap::iterator>(getAtomCount(), "molecule::MinimiseConstrainedPotential: Params.*StepList");
     340  Params.PermutationMap = Calloc<atom*>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.**PermutationMap");
     341  Params.DistanceList = Malloc<DistanceMap*>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.**DistanceList");
     342  Params.DistanceIterators = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*DistanceIterators");
     343  Params.DoubleList = Calloc<int>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*DoubleList");
     344  Params.StepList = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*StepList");
    336345  int round;
    337   atom *Sprinter = NULL;
     346  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
    338347  DistanceMap::iterator Rider, Strider;
    339348
     
    362371    DoLog(2) && (Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl);
    363372    OlderPotential = OldPotential;
    364     molecule::const_iterator iter;
    365373    do {
    366       iter = begin();
    367       for (; iter != end(); ++iter) {
    368         PrintPermutationMap(getAtomCount(), Params);
    369         Sprinter = Params.DistanceIterators[(*iter)->nr]->second;   // store initial partner
    370         Strider = Params.DistanceIterators[(*iter)->nr];  //remember old iterator
    371         Params.DistanceIterators[(*iter)->nr] = Params.StepList[(*iter)->nr];
    372         if (Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->end()) {// stop, before we run through the list and still on
    373           Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->begin();
     374      Walker = start;
     375      while (Walker->next != end) { // pick one
     376        Walker = Walker->next;
     377        PrintPermutationMap(AtomCount, Params);
     378        Sprinter = Params.DistanceIterators[Walker->nr]->second;   // store initial partner
     379        Strider = Params.DistanceIterators[Walker->nr];  //remember old iterator
     380        Params.DistanceIterators[Walker->nr] = Params.StepList[Walker->nr];
     381        if (Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
     382          Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->begin();
    374383          break;
    375384        }
    376         //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->nr]->second << "." << endl;
     385        //Log() << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
    377386        // find source of the new target
    378         molecule::const_iterator runner = begin();
    379         for (; runner != end(); ++runner) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
    380           if (Params.PermutationMap[(*runner)->nr] == Params.DistanceIterators[(*iter)->nr]->second) {
    381             //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->nr] << "." << endl;
     387        Runner = start->next;
     388        while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
     389          if (Params.PermutationMap[Runner->nr] == Params.DistanceIterators[Walker->nr]->second) {
     390            //Log() << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
    382391            break;
    383392          }
     393          Runner = Runner->next;
    384394        }
    385         if (runner != end()) { // we found the other source
     395        if (Runner != end) { // we found the other source
    386396          // then look in its distance list for Sprinter
    387           Rider = Params.DistanceList[(*runner)->nr]->begin();
    388           for (; Rider != Params.DistanceList[(*runner)->nr]->end(); Rider++)
     397          Rider = Params.DistanceList[Runner->nr]->begin();
     398          for (; Rider != Params.DistanceList[Runner->nr]->end(); Rider++)
    389399            if (Rider->second == Sprinter)
    390400              break;
    391           if (Rider != Params.DistanceList[(*runner)->nr]->end()) { // if we have found one
    392             //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->nr] << "/" << *Rider->second << "." << endl;
     401          if (Rider != Params.DistanceList[Runner->nr]->end()) { // if we have found one
     402            //Log() << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl;
    393403            // exchange both
    394             Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // put next farther distance into PermutationMap
    395             Params.PermutationMap[(*runner)->nr] = Sprinter;  // and hand the old target to its respective owner
    396             PrintPermutationMap(getAtomCount(), Params);
     404            Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
     405            Params.PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
     406            PrintPermutationMap(AtomCount, Params);
    397407            // calculate the new potential
    398408            //Log() << Verbose(2) << "Checking new potential ..." << endl;
     
    400410            if (Potential > OldPotential) { // we made everything worse! Undo ...
    401411              //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
    402               //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *Params.DistanceIterators[(*runner)->nr]->second << "." << endl;
     412              //Log() << Verbose(3) << "Setting " << *Runner << "'s source to " << *Params.DistanceIterators[Runner->nr]->second << "." << endl;
    403413              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
    404               Params.PermutationMap[(*runner)->nr] = Params.DistanceIterators[(*runner)->nr]->second;
     414              Params.PermutationMap[Runner->nr] = Params.DistanceIterators[Runner->nr]->second;
    405415              // Undo for Walker
    406               Params.DistanceIterators[(*iter)->nr] = Strider;  // take next farther distance target
    407               //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *Params.DistanceIterators[(*iter)->nr]->second << "." << endl;
    408               Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second;
     416              Params.DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
     417              //Log() << Verbose(3) << "Setting " << *Walker << "'s source to " << *Params.DistanceIterators[Walker->nr]->second << "." << endl;
     418              Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second;
    409419            } else {
    410               Params.DistanceIterators[(*runner)->nr] = Rider;  // if successful also move the pointer in the iterator list
     420              Params.DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
    411421              DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl);
    412422              OldPotential = Potential;
     
    418428            //Log() << Verbose(0) << endl;
    419429          } else {
    420             DoeLog(1) && (eLog()<< Verbose(1) << **runner << " was not the owner of " << *Sprinter << "!" << endl);
     430            DoeLog(1) && (eLog()<< Verbose(1) << *Runner << " was not the owner of " << *Sprinter << "!" << endl);
    421431            exit(255);
    422432          }
    423433        } else {
    424           Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // new target has no source!
     434          Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // new target has no source!
    425435        }
    426         Params.StepList[(*iter)->nr]++; // take next farther distance target
     436        Params.StepList[Walker->nr]++; // take next farther distance target
    427437      }
    428     } while (++iter != end());
     438    } while (Walker->next != end);
    429439  } while ((OlderPotential - OldPotential) > 1e-3);
    430440  DoLog(1) && (Log() << Verbose(1) << "done." << endl);
     
    432442
    433443  /// free memory and return with evaluated potential
    434   for (int i=getAtomCount(); i--;)
     444  for (int i=AtomCount; i--;)
    435445    Params.DistanceList[i]->clear();
    436446  Free(&Params.DistanceList);
     
    473483  // Get the Permutation Map by MinimiseConstrainedPotential
    474484  atom **PermutationMap = NULL;
    475   atom *Sprinter = NULL;
     485  atom *Walker = NULL, *Sprinter = NULL;
    476486  if (!MapByIdentity)
    477487    MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
    478488  else {
    479     PermutationMap = Malloc<atom *>(getAtomCount(), "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
     489    PermutationMap = Malloc<atom *>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
    480490    SetIndexedArrayForEachAtomTo( PermutationMap, &atom::nr );
    481491  }
     
    492502    mol = World::getInstance().createMolecule();
    493503    MoleculePerStep->insert(mol);
    494     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     504    Walker = start;
     505    while (Walker->next != end) {
     506      Walker = Walker->next;
    495507      // add to molecule list
    496       Sprinter = mol->AddCopyAtom((*iter));
     508      Sprinter = mol->AddCopyAtom(Walker);
    497509      for (int n=NDIM;n--;) {
    498         Sprinter->x[n] = (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
     510        Sprinter->x[n] = Walker->Trajectory.R.at(startstep)[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep)[n] - Walker->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
    499511        // add to Trajectories
    500512        //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    501513        if (step < MaxSteps) {
    502           (*iter)->Trajectory.R.at(step)[n] = (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
    503           (*iter)->Trajectory.U.at(step)[n] = 0.;
    504           (*iter)->Trajectory.F.at(step)[n] = 0.;
     514          Walker->Trajectory.R.at(step)[n] = Walker->Trajectory.R.at(startstep)[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep)[n] - Walker->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
     515          Walker->Trajectory.U.at(step)[n] = 0.;
     516          Walker->Trajectory.F.at(step)[n] = 0.;
    505517        }
    506518      }
     
    510522
    511523  // store the list to single step files
    512   int *SortIndex = Malloc<int>(getAtomCount(), "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
    513   for (int i=getAtomCount(); i--; )
     524  int *SortIndex = Malloc<int>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
     525  for (int i=AtomCount; i--; )
    514526    SortIndex[i] = i;
    515527  status = MoleculePerStep->OutputConfigForListOfFragments(&configuration, SortIndex);
     
    555567      return false;
    556568    }
    557     if (Force.RowCounter[0] != getAtomCount()) {
    558       DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
     569    if (Force.RowCounter[0] != AtomCount) {
     570      DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << AtomCount << "." << endl);
    559571      performCriticalExit();
    560572      return false;
     
    562574    // correct Forces
    563575    Velocity.Zero();
    564     for(int i=0;i<getAtomCount();i++)
     576    for(int i=0;i<AtomCount;i++)
    565577      for(int d=0;d<NDIM;d++) {
    566578        Velocity[d] += Force.Matrix[0][i][d+5];
    567579      }
    568     for(int i=0;i<getAtomCount();i++)
     580    for(int i=0;i<AtomCount;i++)
    569581      for(int d=0;d<NDIM;d++) {
    570         Force.Matrix[0][i][d+5] -= Velocity[d]/static_cast<double>(getAtomCount());
     582        Force.Matrix[0][i][d+5] -= Velocity[d]/(double)AtomCount;
    571583      }
    572584    // solve a constrained potential if we are meant to
     
    671683      delta_alpha = 0.;
    672684      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
    673       delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
     685      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
    674686      configuration.alpha += delta_alpha*configuration.Deltat;
    675687      DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl);
  • src/molecule_fragmentation.cpp

    ra7b761b r8f215d  
    3939  int FragmentCount;
    4040  // get maximum bond degree
    41   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    42     c = ((*iter)->ListOfBonds.size() > c) ? (*iter)->ListOfBonds.size() : c;
     41  atom *Walker = start;
     42  while (Walker->next != end) {
     43    Walker = Walker->next;
     44    c = (Walker->ListOfBonds.size() > c) ? Walker->ListOfBonds.size() : c;
    4345  }
    4446  FragmentCount = NoNonHydrogen*(1 << (c*order));
     
    357359map<double, pair<int,int> >  * ReMapAdaptiveCriteriaListToValue(map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol)
    358360{
    359   atom *Walker = NULL;
     361  atom *Walker = mol->start;
    360362  map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ;
    361363  DoLog(1) && (Log() << Verbose(1) << "Root candidate list is: " << endl);
     
    389391bool MarkUpdateCandidates(bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol)
    390392{
    391   atom *Walker = NULL;
     393  atom *Walker = mol->start;
    392394  int No = -1;
    393395  bool status = false;
     
    433435bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path)
    434436{
     437  atom *Walker = start;
    435438  bool status = false;
    436439
    437440  // initialize mask list
    438   for(int i=getAtomCount();i--;)
     441  for(int i=AtomCount;i--;)
    439442    AtomMask[i] = false;
    440443
    441444  if (Order < 0) { // adaptive increase of BondOrder per site
    442     if (AtomMask[getAtomCount()] == true)  // break after one step
     445    if (AtomMask[AtomCount] == true)  // break after one step
    443446      return false;
    444447
     
    454457    if (AdaptiveCriteriaList->empty()) {
    455458      DoeLog(2) && (eLog()<< Verbose(2) << "Unable to parse file, incrementing all." << endl);
    456       for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     459      while (Walker->next != end) {
     460        Walker = Walker->next;
    457461    #ifdef ADDHYDROGEN
    458         if ((*iter)->type->Z != 1) // skip hydrogen
     462        if (Walker->type->Z != 1) // skip hydrogen
    459463    #endif
    460464        {
    461           AtomMask[(*iter)->nr] = true;  // include all (non-hydrogen) atoms
     465          AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
    462466          status = true;
    463467        }
     
    474478    Free(&FinalRootCandidates);
    475479  } else { // global increase of Bond Order
    476     for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     480    while (Walker->next != end) {
     481      Walker = Walker->next;
    477482  #ifdef ADDHYDROGEN
    478       if ((*iter)->type->Z != 1) // skip hydrogen
     483      if (Walker->type->Z != 1) // skip hydrogen
    479484  #endif
    480485      {
    481         AtomMask[(*iter)->nr] = true;  // include all (non-hydrogen) atoms
    482         if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->nr]))
     486        AtomMask[Walker->nr] = true;  // include all (non-hydrogen) atoms
     487        if ((Order != 0) && (Walker->AdaptiveOrder < Order)) // && (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]))
    483488          status = true;
    484489      }
    485490    }
    486     if ((!Order) && (!AtomMask[getAtomCount()]))  // single stepping, just check
     491    if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
    487492      status = true;
    488493
     
    495500  }
    496501
    497   PrintAtomMask(AtomMask, getAtomCount()); // for debugging
     502  PrintAtomMask(AtomMask, AtomCount); // for debugging
    498503
    499504  return status;
     
    511516    return false;
    512517  }
    513   SortIndex = Malloc<int>(getAtomCount(), "molecule::CreateMappingLabelsToConfigSequence: *SortIndex");
    514   for(int i=getAtomCount();i--;)
     518  SortIndex = Malloc<int>(AtomCount, "molecule::CreateMappingLabelsToConfigSequence: *SortIndex");
     519  for(int i=AtomCount;i--;)
    515520    SortIndex[i] = -1;
    516521
     
    519524
    520525  return true;
    521 };
    522 
    523 
    524 
    525 /** Creates a lookup table for true father's Atom::Nr -> atom ptr.
    526  * \param *start begin of list (STL iterator, i.e. first item)
    527  * \paran *end end of list (STL iterator, i.e. one past last item)
    528  * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
    529  * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
    530  * \return true - success, false - failure
    531  */
    532 bool molecule::CreateFatherLookupTable(atom **&LookupTable, int count)
    533 {
    534   bool status = true;
    535   int AtomNo;
    536 
    537   if (LookupTable != NULL) {
    538     Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl;
    539     return false;
    540   }
    541 
    542   // count them
    543   if (count == 0) {
    544     for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
    545       count = (count < (*iter)->GetTrueFather()->nr) ? (*iter)->GetTrueFather()->nr : count;
    546     }
    547   }
    548   if (count <= 0) {
    549     Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl;
    550     return false;
    551   }
    552 
    553   // allocate and fill
    554   LookupTable = Calloc<atom *>(count, "CreateFatherLookupTable - **LookupTable");
    555   if (LookupTable == NULL) {
    556     eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl;
    557     performCriticalExit();
    558     status = false;
    559   } else {
    560     for (molecule::iterator iter = begin(); iter != end(); ++iter) {
    561       AtomNo = (*iter)->GetTrueFather()->nr;
    562       if ((AtomNo >= 0) && (AtomNo < count)) {
    563         //*out << "Setting LookupTable[" << AtomNo << "] to " << *(*iter) << endl;
    564         LookupTable[AtomNo] = (*iter);
    565       } else {
    566         Log() << Verbose(0) << "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << ")." << endl;
    567         status = false;
    568         break;
    569       }
    570     }
    571   }
    572 
    573   return status;
    574526};
    575527
     
    596548  MoleculeListClass *BondFragments = NULL;
    597549  int *SortIndex = NULL;
    598   int *MinimumRingSize = new int[getAtomCount()];
     550  int *MinimumRingSize = new int[AtomCount];
    599551  int FragmentCounter;
    600552  MoleculeLeafClass *MolecularWalker = NULL;
     
    624576
    625577  // create lookup table for Atom::nr
    626   FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(ListOfAtoms, getAtomCount());
     578  FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(start, end, ListOfAtoms, AtomCount);
    627579
    628580  // === compare it with adjacency file ===
     
    634586
    635587  // analysis of the cycles (print rings, get minimum cycle length) for each subgraph
    636   for(int i=getAtomCount();i--;)
    637     MinimumRingSize[i] = getAtomCount();
     588  for(int i=AtomCount;i--;)
     589    MinimumRingSize[i] = AtomCount;
    638590  MolecularWalker = Subgraphs;
    639591  FragmentCounter = 0;
     
    672624  // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
    673625  KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
    674   AtomMask = new bool[getAtomCount()+1];
    675   AtomMask[getAtomCount()] = false;
     626  AtomMask = new bool[AtomCount+1];
     627  AtomMask[AtomCount] = false;
    676628  FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
    677629  while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) {
    678630    FragmentationToDo = FragmentationToDo || CheckOrder;
    679     AtomMask[getAtomCount()] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
     631    AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
    680632    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
    681633    Subgraphs->next->FillRootStackForSubgraphs(RootStack, AtomMask, (FragmentCounter = 0));
     
    816768bool molecule::ParseOrderAtSiteFromFile(char *path)
    817769{
    818   unsigned char *OrderArray = Calloc<unsigned char>(getAtomCount(), "molecule::ParseOrderAtSiteFromFile - *OrderArray");
    819   bool *MaxArray = Calloc<bool>(getAtomCount(), "molecule::ParseOrderAtSiteFromFile - *MaxArray");
     770  unsigned char *OrderArray = Calloc<unsigned char>(AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
     771  bool *MaxArray = Calloc<bool>(AtomCount, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
    820772  bool status;
    821773  int AtomNr, value;
     
    920872  atom *OtherFather = NULL;
    921873  atom *FatherOfRunner = NULL;
    922 
    923 #ifdef ADDHYDROGEN
    924   molecule::const_iterator runner;
    925 #endif
    926   // we increment the iter just before skipping the hydrogen
    927   for (molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end();) {
     874  Leaf->CountAtoms();
     875
     876  atom *Runner = Leaf->start;
     877  while (Runner->next != Leaf->end) {
     878    Runner = Runner->next;
    928879    LonelyFlag = true;
    929     FatherOfRunner = (*iter)->father;
    930     ASSERT(FatherOfRunner,"Atom without father found");
     880    FatherOfRunner = Runner->father;
    931881    if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
    932882      // create all bonds
     
    939889//            Log() << Verbose(3) << "Adding Bond: ";
    940890//            Log() << Verbose(0) <<
    941             Leaf->AddBond((*iter), SonList[OtherFather->nr], (*BondRunner)->BondDegree);
     891            Leaf->AddBond(Runner, SonList[OtherFather->nr], (*BondRunner)->BondDegree);
    942892//            Log() << Verbose(0) << "." << endl;
    943             //NumBonds[(*iter)->nr]++;
     893            //NumBonds[Runner->nr]++;
    944894          } else {
    945895//            Log() << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
     
    949899//          Log() << Verbose(0) << ", who has no son in this fragment molecule." << endl;
    950900#ifdef ADDHYDROGEN
    951           //Log() << Verbose(3) << "Adding Hydrogen to " << (*iter)->Name << " and a bond in between." << endl;
    952           if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
     901          //Log() << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
     902          if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), Runner, FatherOfRunner, OtherFather, IsAngstroem))
    953903            exit(1);
    954904#endif
    955           //NumBonds[(*iter)->nr] += Binder->BondDegree;
     905          //NumBonds[Runner->nr] += Binder->BondDegree;
    956906        }
    957907      }
    958908    } else {
    959     DoeLog(1) && (eLog()<< Verbose(1) << "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl);
    960     }
    961     if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
    962       DoLog(0) && (Log() << Verbose(0) << **iter << "has got bonds only to hydrogens!" << endl);
    963     }
    964     ++iter;
     909      DoeLog(1) && (eLog()<< Verbose(1) << "Son " << Runner->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl);
     910    }
     911    if ((LonelyFlag) && (Leaf->AtomCount > 1)) {
     912      DoLog(0) && (Log() << Verbose(0) << *Runner << "has got bonds only to hydrogens!" << endl);
     913    }
    965914#ifdef ADDHYDROGEN
    966     while ((iter != Leaf->end()) && ((*iter)->type->Z == 1)){ // skip added hydrogen
    967       iter++;
    968     }
     915    while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
     916      Runner = Runner->next;
    969917#endif
    970918  }
     
    981929molecule * molecule::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
    982930{
    983   atom **SonList = Calloc<atom*>(getAtomCount(), "molecule::StoreFragmentFromStack: **SonList");
     931  atom **SonList = Calloc<atom*>(AtomCount, "molecule::StoreFragmentFromStack: **SonList");
    984932  molecule *Leaf = World::getInstance().createMolecule();
    985933
     
    16081556  FragmentSearch.FragmentSet = new KeySet;
    16091557  FragmentSearch.Root = FindAtom(RootKeyNr);
    1610   FragmentSearch.ShortestPathList = Malloc<int>(getAtomCount(), "molecule::PowerSetGenerator: *ShortestPathList");
    1611   for (int i=getAtomCount();i--;) {
     1558  FragmentSearch.ShortestPathList = Malloc<int>(AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
     1559  for (int i=AtomCount;i--;) {
    16121560    FragmentSearch.ShortestPathList[i] = -1;
    16131561  }
    16141562
    16151563  // Construct the complete KeySet which we need for topmost level only (but for all Roots)
     1564  atom *Walker = start;
    16161565  KeySet CompleteMolecule;
    1617   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    1618     CompleteMolecule.insert((*iter)->GetTrueFather()->nr);
     1566  while (Walker->next != end) {
     1567    Walker = Walker->next;
     1568    CompleteMolecule.insert(Walker->GetTrueFather()->nr);
    16191569  }
    16201570
     
    16271577    RootKeyNr = RootStack.front();
    16281578    RootStack.pop_front();
    1629     atom *Walker = FindAtom(RootKeyNr);
     1579    Walker = FindAtom(RootKeyNr);
    16301580    // check cyclic lengths
    16311581    //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
     
    17141664  Vector Translationvector;
    17151665  //class StackClass<atom *> *CompStack = NULL;
    1716   class StackClass<atom *> *AtomStack = new StackClass<atom *>(getAtomCount());
     1666  class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
    17171667  bool flag = true;
    17181668
    17191669  DoLog(2) && (Log() << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl);
    17201670
    1721   ColorList = Calloc<enum Shading>(getAtomCount(), "molecule::ScanForPeriodicCorrection: *ColorList");
     1671  ColorList = Calloc<enum Shading>(AtomCount, "molecule::ScanForPeriodicCorrection: *ColorList");
    17221672  while (flag) {
    17231673    // remove bonds that are beyond bonddistance
     
    17511701      Log() << Verbose(0) << Translationvector <<  endl;
    17521702      // apply to all atoms of first component via BFS
    1753       for (int i=getAtomCount();i--;)
     1703      for (int i=AtomCount;i--;)
    17541704        ColorList[i] = white;
    17551705      AtomStack->Push(Binder->leftatom);
  • src/molecule_geometry.cpp

    ra7b761b r8f215d  
    6969
    7070//  Log() << Verbose(3) << "Begin of CenterEdge." << endl;
    71   molecule::const_iterator iter = begin();  // start at first in list
    72   if (iter != end()) { //list not empty?
     71  atom *ptr = start->next;  // start at first in list
     72  if (ptr != end) {  //list not empty?
    7373    for (int i=NDIM;i--;) {
    74       max->at(i) = (*iter)->x[i];
    75       min->at(i) = (*iter)->x[i];
    76     }
    77     for (; iter != end(); ++iter) {// continue with second if present
    78       //(*iter)->Output(1,1,out);
     74      max->at(i) = ptr->x[i];
     75      min->at(i) = ptr->x[i];
     76    }
     77    while (ptr->next != end) {  // continue with second if present
     78      ptr = ptr->next;
     79      //ptr->Output(1,1,out);
    7980      for (int i=NDIM;i--;) {
    80         max->at(i) = (max->at(i) < (*iter)->x[i]) ? (*iter)->x[i] : max->at(i);
    81         min->at(i) = (min->at(i) > (*iter)->x[i]) ? (*iter)->x[i] : min->at(i);
     81        max->at(i) = (max->at(i) < ptr->x[i]) ? ptr->x[i] : max->at(i);
     82        min->at(i) = (min->at(i) > ptr->x[i]) ? ptr->x[i] : min->at(i);
    8283      }
    8384    }
     
    103104{
    104105  int Num = 0;
    105   molecule::const_iterator iter = begin();  // start at first in list
     106  atom *ptr = start;  // start at first in list
    106107
    107108  Center.Zero();
    108109
    109   if (iter != end()) {   //list not empty?
    110     for (; iter != end(); ++iter) {  // continue with second if present
     110  if (ptr->next != end) {   //list not empty?
     111    while (ptr->next != end) {  // continue with second if present
     112      ptr = ptr->next;
    111113      Num++;
    112       Center += (*iter)->x;
     114      Center += ptr->x;
    113115    }
    114116    Center.Scale(-1./Num); // divide through total number (and sign for direction)
     
    123125Vector * molecule::DetermineCenterOfAll() const
    124126{
    125   molecule::const_iterator iter = begin();  // start at first in list
     127  atom *ptr = start->next;  // start at first in list
    126128  Vector *a = new Vector();
    127129  Vector tmp;
     
    130132  a->Zero();
    131133
    132   if (iter != end()) {   //list not empty?
    133     for (; iter != end(); ++iter) {  // continue with second if present
     134  if (ptr != end) {   //list not empty?
     135    while (ptr->next != end) {  // continue with second if present
     136      ptr = ptr->next;
    134137      Num += 1.;
    135       tmp = (*iter)->x;
     138      tmp = ptr->x;
    136139      (*a) += tmp;
    137140    }
     
    147150Vector * molecule::DetermineCenterOfGravity()
    148151{
    149   molecule::const_iterator iter = begin();  // start at first in list
     152  atom *ptr = start->next;  // start at first in list
    150153  Vector *a = new Vector();
    151154  Vector tmp;
     
    154157  a->Zero();
    155158
    156   if (iter != end()) {   //list not empty?
    157     for (; iter != end(); ++iter) {  // continue with second if present
    158       Num += (*iter)->type->mass;
    159       tmp = (*iter)->type->mass * (*iter)->x;
     159  if (ptr != end) {   //list not empty?
     160    while (ptr->next != end) {  // continue with second if present
     161      ptr = ptr->next;
     162      Num += ptr->type->mass;
     163      tmp = ptr->type->mass * ptr->x;
    160164      (*a) += tmp;
    161165    }
     
    198202void molecule::Scale(const double ** const factor)
    199203{
    200   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     204  atom *ptr = start;
     205
     206  while (ptr->next != end) {
     207    ptr = ptr->next;
    201208    for (int j=0;j<MDSteps;j++)
    202       (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
    203     (*iter)->x.ScaleAll(*factor);
     209      ptr->Trajectory.R.at(j).ScaleAll(*factor);
     210    ptr->x.ScaleAll(*factor);
    204211  }
    205212};
     
    210217void molecule::Translate(const Vector *trans)
    211218{
    212   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     219  atom *ptr = start;
     220
     221  while (ptr->next != end) {
     222    ptr = ptr->next;
    213223    for (int j=0;j<MDSteps;j++)
    214       (*iter)->Trajectory.R.at(j) += (*trans);
    215     (*iter)->x += (*trans);
     224      ptr->Trajectory.R.at(j) += (*trans);
     225    ptr->x += (*trans);
    216226  }
    217227};
     
    249259void molecule::DeterminePeriodicCenter(Vector &center)
    250260{
     261  atom *Walker = start;
    251262  double * const cell_size = World::getInstance().getDomain();
    252263  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
     
    259270    Center.Zero();
    260271    flag = true;
    261     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     272    while (Walker->next != end) {
     273      Walker = Walker->next;
    262274#ifdef ADDHYDROGEN
    263       if ((*iter)->type->Z != 1) {
     275      if (Walker->type->Z != 1) {
    264276#endif
    265         Testvector = (*iter)->x;
     277        Testvector = Walker->x;
    266278        Testvector.MatrixMultiplication(inversematrix);
    267279        Translationvector.Zero();
    268         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    269          if ((*iter)->nr < (*Runner)->GetOtherAtom((*iter))->nr) // otherwise we shift one to, the other fro and gain nothing
     280        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
     281         if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
    270282            for (int j=0;j<NDIM;j++) {
    271               tmp = (*iter)->x[j] - (*Runner)->GetOtherAtom(*iter)->x[j];
     283              tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
    272284              if ((fabs(tmp)) > BondDistance) {
    273285                flag = false;
    274                 DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << (*iter)->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
     286                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
    275287                if (tmp > 0)
    276288                  Translationvector[j] -= 1.;
     
    286298#ifdef ADDHYDROGEN
    287299        // now also change all hydrogens
    288         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    289           if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
    290             Testvector = (*Runner)->GetOtherAtom((*iter))->x;
     300        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
     301          if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
     302            Testvector = (*Runner)->GetOtherAtom(Walker)->x;
    291303            Testvector.MatrixMultiplication(inversematrix);
    292304            Testvector += Translationvector;
     
    303315  Free(&inversematrix);
    304316
    305   Center.Scale(1./static_cast<double>(getAtomCount()));
     317  Center.Scale(1./(double)AtomCount);
    306318};
    307319
     
    313325void molecule::PrincipalAxisSystem(bool DoRotate)
    314326{
     327  atom *ptr = start;  // start at first in list
    315328  double InertiaTensor[NDIM*NDIM];
    316329  Vector *CenterOfGravity = DetermineCenterOfGravity();
     
    323336
    324337  // sum up inertia tensor
    325   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    326     Vector x = (*iter)->x;
     338  while (ptr->next != end) {
     339    ptr = ptr->next;
     340    Vector x = ptr->x;
    327341    //x.SubtractVector(CenterOfGravity);
    328     InertiaTensor[0] += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
    329     InertiaTensor[1] += (*iter)->type->mass*(-x[0]*x[1]);
    330     InertiaTensor[2] += (*iter)->type->mass*(-x[0]*x[2]);
    331     InertiaTensor[3] += (*iter)->type->mass*(-x[1]*x[0]);
    332     InertiaTensor[4] += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
    333     InertiaTensor[5] += (*iter)->type->mass*(-x[1]*x[2]);
    334     InertiaTensor[6] += (*iter)->type->mass*(-x[2]*x[0]);
    335     InertiaTensor[7] += (*iter)->type->mass*(-x[2]*x[1]);
    336     InertiaTensor[8] += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
     342    InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
     343    InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
     344    InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
     345    InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
     346    InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
     347    InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
     348    InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
     349    InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
     350    InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
    337351  }
    338352  // print InertiaTensor for debugging
     
    372386
    373387    // sum up inertia tensor
    374     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    375       Vector x = (*iter)->x;
    376       InertiaTensor[0] += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
    377       InertiaTensor[1] += (*iter)->type->mass*(-x[0]*x[1]);
    378       InertiaTensor[2] += (*iter)->type->mass*(-x[0]*x[2]);
    379       InertiaTensor[3] += (*iter)->type->mass*(-x[1]*x[0]);
    380       InertiaTensor[4] += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
    381       InertiaTensor[5] += (*iter)->type->mass*(-x[1]*x[2]);
    382       InertiaTensor[6] += (*iter)->type->mass*(-x[2]*x[0]);
    383       InertiaTensor[7] += (*iter)->type->mass*(-x[2]*x[1]);
    384       InertiaTensor[8] += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
     388    ptr = start;
     389    while (ptr->next != end) {
     390      ptr = ptr->next;
     391      Vector x = ptr->x;
     392      //x.SubtractVector(CenterOfGravity);
     393      InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
     394      InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
     395      InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
     396      InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
     397      InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
     398      InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
     399      InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
     400      InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
     401      InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
    385402    }
    386403    // print InertiaTensor for debugging
     
    406423void molecule::Align(Vector *n)
    407424{
     425  atom *ptr = start;
    408426  double alpha, tmp;
    409427  Vector z_axis;
     
    416434  alpha = atan(-n->at(0)/n->at(2));
    417435  DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
    418   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    419     tmp = (*iter)->x[0];
    420     (*iter)->x[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x[2];
    421     (*iter)->x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x[2];
     436  while (ptr->next != end) {
     437    ptr = ptr->next;
     438    tmp = ptr->x[0];
     439    ptr->x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
     440    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
    422441    for (int j=0;j<MDSteps;j++) {
    423       tmp = (*iter)->Trajectory.R.at(j)[0];
    424       (*iter)->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
    425       (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
     442      tmp = ptr->Trajectory.R.at(j)[0];
     443      ptr->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
     444      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
    426445    }
    427446  }
     
    433452
    434453  // rotate on z-y plane
     454  ptr = start;
    435455  alpha = atan(-n->at(1)/n->at(2));
    436456  DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
    437   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    438     tmp = (*iter)->x[1];
    439     (*iter)->x[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x[2];
    440     (*iter)->x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x[2];
     457  while (ptr->next != end) {
     458    ptr = ptr->next;
     459    tmp = ptr->x[1];
     460    ptr->x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
     461    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
    441462    for (int j=0;j<MDSteps;j++) {
    442       tmp = (*iter)->Trajectory.R.at(j)[1];
    443       (*iter)->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
    444       (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
     463      tmp = ptr->Trajectory.R.at(j)[1];
     464      ptr->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
     465      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
    445466    }
    446467  }
     
    466487  Vector a,b,c,d;
    467488  struct lsq_params *par = (struct lsq_params *)params;
     489  atom *ptr = par->mol->start;
    468490
    469491  // initialize vectors
     
    475497  b[2] = gsl_vector_get(x,5);
    476498  // go through all atoms
    477   for (molecule::const_iterator iter = par->mol->begin(); iter != par->mol->end(); ++iter) {
    478     if ((*iter)->type == ((struct lsq_params *)params)->type) { // for specific type
    479       c = (*iter)->x - a;
     499  while (ptr != par->mol->end) {
     500    ptr = ptr->next;
     501    if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
     502      c = ptr->x - a;
    480503      t = c.ScalarProduct(b);           // get direction parameter
    481504      d = t*b;       // and create vector
  • src/molecule_graph.cpp

    ra7b761b r8f215d  
    1919#include "World.hpp"
    2020#include "Helpers/fast_functions.hpp"
    21 #include "Helpers/Assert.hpp"
    22 
    2321
    2422struct BFSAccounting
     
    7674      flip(atom1, atom2);
    7775    Walker = FindAtom(atom1);
    78     ASSERT(Walker,"Could not find an atom with the ID given in dbond file");
    7976    OtherWalker = FindAtom(atom2);
    80     ASSERT(OtherWalker,"Could not find an atom with the ID given in dbond file");
    8177    AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
    8278  }
     
    107103  atom *Walker = NULL;
    108104  atom *OtherWalker = NULL;
     105  atom **AtomMap = NULL;
    109106  int n[NDIM];
    110107  double MinDistance, MaxDistance;
     
    131128
    132129  // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
    133   DoLog(1) && (Log() << Verbose(1) << "AtomCount " << getAtomCount() << " and bonddistance is " << bonddistance << "." << endl);
    134 
    135   if ((getAtomCount() > 1) && (bonddistance > 1.)) {
     130  CountAtoms();
     131  DoLog(1) && (Log() << Verbose(1) << "AtomCount " << AtomCount << " and bonddistance is " << bonddistance << "." << endl);
     132
     133  if ((AtomCount > 1) && (bonddistance > 1.)) {
    136134    DoLog(2) && (Log() << Verbose(2) << "Creating Linked Cell structure ... " << endl);
    137135    LC = new LinkedCell(this, bonddistance);
     
    139137    // create a list to map Tesselpoint::nr to atom *
    140138    DoLog(2) && (Log() << Verbose(2) << "Creating TesselPoint to atom map ... " << endl);
    141 
    142     // set numbers for atoms that can later be used
    143     int i=0;
    144     for(internal_iterator iter = atoms.begin();iter!= atoms.end(); ++iter){
    145       (*iter)->nr = i++;
     139    AtomMap = Calloc<atom *> (AtomCount, "molecule::CreateAdjacencyList - **AtomCount");
     140    Walker = start;
     141    while (Walker->next != end) {
     142      Walker = Walker->next;
     143      AtomMap[Walker->nr] = Walker;
    146144    }
    147145
     
    155153          if (List != NULL) {
    156154            for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    157               Walker = dynamic_cast<atom*>(*Runner);
    158               ASSERT(Walker,"Tesselpoint that was not an atom retrieved from LinkedNode");
    159               //Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
     155              Walker = AtomMap[(*Runner)->nr];
     156//              Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
    160157              // 3c. check for possible bond between each atom in this and every one in the 27 cells
    161158              for (n[0] = -1; n[0] <= 1; n[0]++)
     
    167164                      for (LinkedCell::LinkedNodes::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
    168165                        if ((*OtherRunner)->nr > Walker->nr) {
    169                           OtherWalker = dynamic_cast<atom*>(*OtherRunner);
    170                           ASSERT(OtherWalker,"TesselPoint that was not an atom retrieved from LinkedNode");
    171                           //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
     166                          OtherWalker = AtomMap[(*OtherRunner)->nr];
     167//                          Log() << Verbose(0) << "Current other Atom is " << *OtherWalker << "." << endl;
     168                          const double distance = OtherWalker->x.PeriodicDistanceSquared(Walker->x, cell_size);
     169//                          Log() << Verbose(1) << "Checking distance " << distance << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
    172170                          (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
    173                           const double distance = OtherWalker->x.PeriodicDistanceSquared(Walker->x,cell_size);
    174171                          const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
    175172//                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
     
    191188          }
    192189        }
     190    Free(&AtomMap);
    193191    delete (LC);
    194192    DoLog(1) && (Log() << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl);
     
    201199    ActOnAllAtoms( &atom::OutputBondOfAtom );
    202200  } else
    203     DoLog(1) && (Log() << Verbose(1) << "AtomCount is " << getAtomCount() << ", thus no bonds, no connections!." << endl);
     201    DoLog(1) && (Log() << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl);
    204202  DoLog(0) && (Log() << Verbose(0) << "End of CreateAdjacencyList." << endl);
    205203  if (free_BG)
     
    243241    DoLog(0) && (Log() << Verbose(0) << " done." << endl);
    244242  } else {
    245     DoLog(1) && (Log() << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << getAtomCount() << " atoms." << endl);
     243    DoLog(1) && (Log() << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl);
    246244  }
    247245  DoLog(0) && (Log() << Verbose(0) << No << " bonds could not be corrected." << endl);
     
    463461void DepthFirstSearchAnalysis_Init(struct DFSAccounting &DFS, const molecule * const mol)
    464462{
    465   DFS.AtomStack = new StackClass<atom *> (mol->getAtomCount());
     463  DFS.AtomStack = new StackClass<atom *> (mol->AtomCount);
    466464  DFS.CurrentGraphNr = 0;
    467465  DFS.ComponentNumber = 0;
     
    504502  bond *Binder = NULL;
    505503
    506   if (getAtomCount() == 0)
     504  if (AtomCount == 0)
    507505    return SubGraphs;
    508506  DoLog(0) && (Log() << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl);
    509507  DepthFirstSearchAnalysis_Init(DFS, this);
    510508
    511   for (molecule::const_iterator iter = begin(); iter != end();) {
    512     DFS.Root = *iter;
     509  DFS.Root = start->next;
     510  while (DFS.Root != end) { // if there any atoms at all
    513511    // (1) mark all edges unused, empty stack, set atom->GraphNr = -1 for all
    514512    DFS.AtomStack->ClearStack();
     
    550548
    551549    // step on to next root
    552     while ((iter != end()) && ((*iter)->GraphNr != -1)) {
    553       //Log() << Verbose(1) << "Current next subgraph root candidate is " << (*iter)->Name << "." << endl;
    554       if ((*iter)->GraphNr != -1) // if already discovered, step on
    555         iter++;
     550    while ((DFS.Root != end) && (DFS.Root->GraphNr != -1)) {
     551      //Log() << Verbose(1) << "Current next subgraph root candidate is " << Root->Name << "." << endl;
     552      if (DFS.Root->GraphNr != -1) // if already discovered, step on
     553        DFS.Root = DFS.Root->next;
    556554    }
    557555  }
     
    856854  if (MinRingSize != -1) { // if rings are present
    857855    // go over all atoms
    858     for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    859       Root = *iter;
    860 
    861       if (MinimumRingSize[Root->GetTrueFather()->nr] == mol->getAtomCount()) { // check whether MinimumRingSize is set, if not BFS to next where it is
     856    Root = mol->start;
     857    while (Root->next != mol->end) {
     858      Root = Root->next;
     859
     860      if (MinimumRingSize[Root->GetTrueFather()->nr] == mol->AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
    862861        Walker = Root;
    863862
    864863        //Log() << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
    865         CyclicStructureAnalysis_BFSToNextCycle(Root, Walker, MinimumRingSize, mol->getAtomCount());
     864        CyclicStructureAnalysis_BFSToNextCycle(Root, Walker, MinimumRingSize, mol->AtomCount);
    866865
    867866      }
     
    893892  int MinRingSize = -1;
    894893
    895   InitializeBFSAccounting(BFS, getAtomCount());
     894  InitializeBFSAccounting(BFS, AtomCount);
    896895
    897896  //Log() << Verbose(1) << "Back edge list - ";
     
    11351134    CurrentBondsOfAtom = -1; // we count one too far due to line end
    11361135    // parse into structure
    1137     if ((AtomNr >= 0) && (AtomNr < getAtomCount())) {
     1136    if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
    11381137      Walker = ListOfAtoms[AtomNr];
    11391138      while (!line.eof())
     
    13141313    AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root);
    13151314
    1316   BreadthFirstSearchAdd_Init(BFS, Root, BondOrder, getAtomCount(), AddedAtomList);
     1315  BreadthFirstSearchAdd_Init(BFS, Root, BondOrder, AtomCount, AddedAtomList);
    13171316
    13181317  // and go on ... Queue always contains all lightgray vertices
     
    13701369  // fill parent list with sons
    13711370  DoLog(3) && (Log() << Verbose(3) << "Filling Parent List." << endl);
    1372   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    1373     ParentList[(*iter)->father->nr] = (*iter);
     1371  atom *Walker = mol->start;
     1372  while (Walker->next != mol->end) {
     1373    Walker = Walker->next;
     1374    ParentList[Walker->father->nr] = Walker;
    13741375    // Outputting List for debugging
    1375     DoLog(4) && (Log() << Verbose(4) << "Son[" << (*iter)->father->nr << "] of " << (*iter)->father << " is " << ParentList[(*iter)->father->nr] << "." << endl);
    1376   }
    1377 };
     1376    DoLog(4) && (Log() << Verbose(4) << "Son[" << Walker->father->nr << "] of " << Walker->father << " is " << ParentList[Walker->father->nr] << "." << endl);
     1377  }
     1378
     1379}
     1380;
    13781381
    13791382void BuildInducedSubgraph_Finalize(atom **&ParentList)
     
    13861389{
    13871390  bool status = true;
     1391  atom *Walker = NULL;
    13881392  atom *OtherAtom = NULL;
    13891393  // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
    13901394  DoLog(3) && (Log() << Verbose(3) << "Creating bonds." << endl);
    1391   for (molecule::const_iterator iter = Father->begin(); iter != Father->end(); ++iter) {
    1392     if (ParentList[(*iter)->nr] != NULL) {
    1393       if (ParentList[(*iter)->nr]->father != (*iter)) {
     1395  Walker = Father->start;
     1396  while (Walker->next != Father->end) {
     1397    Walker = Walker->next;
     1398    if (ParentList[Walker->nr] != NULL) {
     1399      if (ParentList[Walker->nr]->father != Walker) {
    13941400        status = false;
    13951401      } else {
    1396         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    1397           OtherAtom = (*Runner)->GetOtherAtom((*iter));
     1402        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
     1403          OtherAtom = (*Runner)->GetOtherAtom(Walker);
    13981404          if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
    1399             DoLog(4) && (Log() << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[(*iter)->nr]->getName() << " and " << ParentList[OtherAtom->nr]->getName() << "." << endl);
    1400             mol->AddBond(ParentList[(*iter)->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
     1405            DoLog(4) && (Log() << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[Walker->nr]->getName() << " and " << ParentList[OtherAtom->nr]->getName() << "." << endl);
     1406            mol->AddBond(ParentList[Walker->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
    14011407          }
    14021408        }
     
    14211427  bool status = true;
    14221428  atom **ParentList = NULL;
     1429
    14231430  DoLog(2) && (Log() << Verbose(2) << "Begin of BuildInducedSubgraph." << endl);
    1424   BuildInducedSubgraph_Init(ParentList, Father->getAtomCount());
     1431  BuildInducedSubgraph_Init(ParentList, Father->AtomCount);
    14251432  BuildInducedSubgraph_FillParentList(this, Father, ParentList);
    14261433  status = BuildInducedSubgraph_CreateBondsFromParent(this, Father, ParentList);
  • src/molecule_pointcloud.cpp

    ra7b761b r8f215d  
    88#include "atom.hpp"
    99#include "config.hpp"
    10 #include "info.hpp"
    1110#include "memoryallocator.hpp"
    1211#include "molecule.hpp"
     
    3231};
    3332
    34 
    35 /** PointCloud implementation of GoPoint
    36  * Uses atoms and STL stuff.
     33/** Return current atom in the list.
     34 * \return pointer to atom or NULL if none present
    3735 */
    38 TesselPoint* molecule::GetPoint() const
     36TesselPoint *molecule::GetPoint() const
    3937{
    40   Info FunctionInfo(__func__);
    41   return (*InternalPointer);
     38  if ((InternalPointer != start) && (InternalPointer != end))
     39    return InternalPointer;
     40  else
     41    return NULL;
    4242};
    4343
    44 /** PointCloud implementation of GoToNext.
    45  * Uses atoms and STL stuff.
     44/** Return pointer to one after last atom in the list.
     45 * \return pointer to end marker
     46 */
     47TesselPoint *molecule::GetTerminalPoint() const
     48{
     49  return end;
     50};
     51
     52/** Return the greatest index of all atoms in the list.
     53 * \return greatest index
     54 */
     55int molecule::GetMaxId() const
     56{
     57  return last_atom;
     58};
     59
     60/** Go to next atom.
     61 * Stops at last one.
    4662 */
    4763void molecule::GoToNext() const
    4864{
    49   Info FunctionInfo(__func__);
    50   if (InternalPointer != atoms.end())
    51     InternalPointer++;
     65  if (InternalPointer != end)
     66    InternalPointer = InternalPointer->next;
    5267};
    5368
    54 /** PointCloud implementation of GoToFirst.
    55  * Uses atoms and STL stuff.
     69/** Go to previous atom.
     70 * Stops at first one.
     71 */
     72void molecule::GoToPrevious() const
     73{
     74  if (InternalPointer->previous != start)
     75    InternalPointer = InternalPointer->previous;
     76};
     77
     78/** Goes to first atom.
    5679 */
    5780void molecule::GoToFirst() const
    5881{
    59   Info FunctionInfo(__func__);
    60   InternalPointer = atoms.begin();
     82  InternalPointer = start->next;
    6183};
    6284
    63 /** PointCloud implementation of IsEmpty.
    64  * Uses atoms and STL stuff.
     85/** Goes to last atom.
     86 */
     87void molecule::GoToLast() const
     88{
     89  InternalPointer = end->previous;
     90};
     91
     92/** Checks whether we have any atoms in molecule.
     93 * \return true - no atoms, false - not empty
    6594 */
    6695bool molecule::IsEmpty() const
    6796{
    68   Info FunctionInfo(__func__);
    69   return (empty());
     97  return (start->next == end);
    7098};
    7199
    72 /** PointCloud implementation of IsLast.
    73  * Uses atoms and STL stuff.
     100/** Checks whether we are at the last atom
     101 * \return true - current atom is last one, false - is not last one
    74102 */
    75103bool molecule::IsEnd() const
    76104{
    77   Info FunctionInfo(__func__);
    78   return (InternalPointer == atoms.end());
     105  return (InternalPointer == end);
    79106};
    80 
    81 int molecule::GetMaxId() const {
    82   return getAtomCount();
    83 }
  • src/molecule_template.hpp

    ra7b761b r8f215d  
    2424template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() ) const
    2525    {
    26   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    27     (((*iter)->node)->*f)();
     26  atom *Walker = start;
     27  while (Walker->next != end) {
     28    Walker = Walker->next;
     29    ((Walker->node)->*f)();
    2830  }
    2931};
    3032template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() const ) const
    3133    {
    32   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    33     (((*iter)->node)->*f)();
     34  atom *Walker = start;
     35  while (Walker->next != end) {
     36    Walker = Walker->next;
     37    ((Walker->node)->*f)();
    3438  }
    3539};
     
    3741template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t ) const
    3842{
    39   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    40     (((*iter)->node)->*f)(t);
     43  atom *Walker = start;
     44  while (Walker->next != end) {
     45    Walker = Walker->next;
     46    ((Walker->node)->*f)(t);
    4147  }
    4248};
    4349template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T) const, T t ) const
    4450{
    45   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    46     (((*iter)->node)->*f)(t);
     51  atom *Walker = start;
     52  while (Walker->next != end) {
     53    Walker = Walker->next;
     54    ((Walker->node)->*f)(t);
    4755  }
    4856};
    4957template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T&), T &t ) const
    5058{
    51   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    52     (((*iter)->node)->*f)(t);
     59  atom *Walker = start;
     60  while (Walker->next != end) {
     61    Walker = Walker->next;
     62    ((Walker->node)->*f)(t);
    5363  }
    5464};
    5565template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T&) const, T &t ) const
    5666{
    57   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    58     (((*iter)->node)->*f)(t);
     67  atom *Walker = start;
     68  while (Walker->next != end) {
     69    Walker = Walker->next;
     70    ((Walker->node)->*f)(t);
    5971  }
    6072};
     
    6274template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const
    6375{
    64   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    65     (((*iter)->node)->*f)(t, u);
     76  atom *Walker = start;
     77  while (Walker->next != end) {
     78    Walker = Walker->next;
     79    ((Walker->node)->*f)(t, u);
    6680  }
    6781};
    6882template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const
    6983{
    70   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    71     (((*iter)->node)->*f)(t, u);
     84  atom *Walker = start;
     85  while (Walker->next != end) {
     86    Walker = Walker->next;
     87    ((Walker->node)->*f)(t, u);
    7288  }
    7389};
     
    7591template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const
    7692{
    77   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    78     (((*iter)->node)->*f)(t, u, v);
     93  atom *Walker = start;
     94  while (Walker->next != end) {
     95    Walker = Walker->next;
     96    ((Walker->node)->*f)(t, u, v);
    7997  }
    8098};
    8199template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V) const, T t, U u, V v) const
    82100{
    83   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    84     (((*iter)->node)->*f)(t, u, v);
     101  atom *Walker = start;
     102  while (Walker->next != end) {
     103    Walker = Walker->next;
     104    ((Walker->node)->*f)(t, u, v);
    85105  }
    86106};
     
    92112{
    93113  res result = 0;
    94   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    95     result += ((*iter)->*f)();
     114  atom *Walker = start;
     115  while (Walker->next != end) {
     116    Walker = Walker->next;
     117    result += (Walker->*f)();
    96118  }
    97119  return result;
     
    100122{
    101123  res result = 0;
    102   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    103     result += ((*iter)->*f)();
     124  atom *Walker = start;
     125  while (Walker->next != end) {
     126    Walker = Walker->next;
     127    result += (Walker->*f)();
    104128  }
    105129  return result;
     
    109133{
    110134  res result = 0;
    111   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    112     result += ((*iter)->*f)(t);
     135  atom *Walker = start;
     136  while (Walker->next != end) {
     137    Walker = Walker->next;
     138    result += (Walker->*f)(t);
    113139  }
    114140  return result;
     
    117143{
    118144  res result = 0;
    119   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    120     result += ((*iter)->*f)(t);
     145  atom *Walker = start;
     146  while (Walker->next != end) {
     147    Walker = Walker->next;
     148    result += (Walker->*f)(t);
    121149  }
    122150  return result;
     
    129157template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *)) const
    130158{
    131   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    132     (*f)((*iter));
     159  atom *Walker = start;
     160  while (Walker->next != end) {
     161    Walker = Walker->next;
     162    (*f)(Walker);
    133163  }
    134164};
    135165template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *) const) const
    136166{
    137   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    138     (*f)((*iter));
     167  atom *Walker = start;
     168  while (Walker->next != end) {
     169    Walker = Walker->next;
     170    (*f)(Walker);
    139171  }
    140172};
     
    145177template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy) const
    146178{
    147   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    148     (copy->*f)((*iter));
     179  atom *Walker = start;
     180  while (Walker->next != end) {
     181    Walker = Walker->next;
     182    (copy->*f)(Walker);
    149183  }
    150184};
    151185template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const, molecule *copy) const
    152186{
    153   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    154     (copy->*f)((*iter));
     187  atom *Walker = start;
     188  while (Walker->next != end) {
     189    Walker = Walker->next;
     190    (copy->*f)(Walker);
    155191  }
    156192};
     
    161197template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () ) const
    162198{
    163   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    164     if (((*iter)->*condition)())
    165       (copy->*f)((*iter));
     199  atom *Walker = start;
     200  while (Walker->next != end) {
     201    Walker = Walker->next;
     202    if ((Walker->*condition)())
     203      (copy->*f)(Walker);
    166204  }
    167205};
    168206template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () const ) const
    169207{
    170   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    171     if (((*iter)->*condition)())
    172       (copy->*f)((*iter));
     208  atom *Walker = start;
     209  while (Walker->next != end) {
     210    Walker = Walker->next;
     211    if ((Walker->*condition)())
     212      (copy->*f)(Walker);
    173213  }
    174214};
    175215template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) () ) const
    176216{
    177   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    178     if (((*iter)->*condition)())
    179       (copy->*f)((*iter));
     217  atom *Walker = start;
     218  while (Walker->next != end) {
     219    Walker = Walker->next;
     220    if ((Walker->*condition)())
     221      (copy->*f)(Walker);
    180222  }
    181223};
    182224template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) () const ) const
    183225{
    184   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    185     if (((*iter)->*condition)())
    186       (copy->*f)((*iter));
     226  atom *Walker = start;
     227  while (Walker->next != end) {
     228    Walker = Walker->next;
     229    if ((Walker->*condition)())
     230      (copy->*f)(Walker);
    187231  }
    188232};
     
    190234template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t ) const
    191235{
    192   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    193     if (((*iter)->*condition)(t))
    194       (copy->*f)((*iter));
     236  atom *Walker = start;
     237  while (Walker->next != end) {
     238    Walker = Walker->next;
     239    if ((Walker->*condition)(t))
     240      (copy->*f)(Walker);
    195241  }
    196242};
    197243template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T) const, T t ) const
    198244{
    199   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    200     if (((*iter)->*condition)(t))
    201       (copy->*f)((*iter));
     245  atom *Walker = start;
     246  while (Walker->next != end) {
     247    Walker = Walker->next;
     248    if ((Walker->*condition)(t))
     249      (copy->*f)(Walker);
    202250  }
    203251};
    204252template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T), T t ) const
    205253{
    206   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    207     if (((*iter)->*condition)(t))
    208       (copy->*f)((*iter));
     254  atom *Walker = start;
     255  while (Walker->next != end) {
     256    Walker = Walker->next;
     257    if ((Walker->*condition)(t))
     258      (copy->*f)(Walker);
    209259  }
    210260};
    211261template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T) const, T t ) const
    212262{
    213   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    214     if (((*iter)->*condition)(t))
    215       (copy->*f)((*iter));
     263  atom *Walker = start;
     264  while (Walker->next != end) {
     265    Walker = Walker->next;
     266    if ((Walker->*condition)(t))
     267      (copy->*f)(Walker);
    216268  }
    217269};
     
    219271template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const
    220272{
    221   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    222     if (((*iter)->*condition)(t,u))
    223       (copy->*f)((*iter));
     273  atom *Walker = start;
     274  while (Walker->next != end) {
     275    Walker = Walker->next;
     276    if ((Walker->*condition)(t,u))
     277      (copy->*f)(Walker);
    224278  }
    225279};
    226280template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const
    227281{
    228   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    229     if (((*iter)->*condition)(t,u))
    230       (copy->*f)((*iter));
     282  atom *Walker = start;
     283  while (Walker->next != end) {
     284    Walker = Walker->next;
     285    if ((Walker->*condition)(t,u))
     286      (copy->*f)(Walker);
    231287  }
    232288};
    233289template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U), T t, U u ) const
    234290{
    235   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    236     if (((*iter)->*condition)(t,u))
    237       (copy->*f)((*iter));
     291  atom *Walker = start;
     292  while (Walker->next != end) {
     293    Walker = Walker->next;
     294    if ((Walker->*condition)(t,u))
     295      (copy->*f)(Walker);
    238296  }
    239297};
    240298template <typename res, typename T, typename U> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U) const, T t, U u ) const
    241299{
    242   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    243     if (((*iter)->*condition)(t,u))
    244       (copy->*f)((*iter));
     300  atom *Walker = start;
     301  while (Walker->next != end) {
     302    Walker = Walker->next;
     303    if ((Walker->*condition)(t,u))
     304      (copy->*f)(Walker);
    245305  }
    246306};
     
    248308template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const
    249309{
    250   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    251     if (((*iter)->*condition)(t,u,v))
    252       (copy->*f)((*iter));
     310  atom *Walker = start;
     311  while (Walker->next != end) {
     312    Walker = Walker->next;
     313    if ((Walker->*condition)(t,u,v))
     314      (copy->*f)(Walker);
    253315  }
    254316};
    255317template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const
    256318{
    257   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    258     if (((*iter)->*condition)(t,u,v))
    259       (copy->*f)((*iter));
     319  atom *Walker = start;
     320  while (Walker->next != end) {
     321    Walker = Walker->next;
     322    if ((Walker->*condition)(t,u,v))
     323      (copy->*f)(Walker);
    260324  }
    261325};
    262326template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U, V), T t, U u, V v ) const
    263327{
    264   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    265     if (((*iter)->*condition)(t,u,v))
    266       (copy->*f)((*iter));
     328  atom *Walker = start;
     329  while (Walker->next != end) {
     330    Walker = Walker->next;
     331    if ((Walker->*condition)(t,u,v))
     332      (copy->*f)(Walker);
    267333  }
    268334};
    269335template <typename res, typename T, typename U, typename V> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T, U, V) const, T t, U u, V v ) const
    270336{
    271   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    272     if (((*iter)->*condition)(t,u,v))
    273       (copy->*f)((*iter));
     337  atom *Walker = start;
     338  while (Walker->next != end) {
     339    Walker = Walker->next;
     340    if ((Walker->*condition)(t,u,v))
     341      (copy->*f)(Walker);
    274342  }
    275343};
     
    280348template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)()) const
    281349{
    282   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    283     ((*iter)->*f)();
     350  atom *Walker = start;
     351  while (Walker->next != end) {
     352    Walker = Walker->next;
     353    (Walker->*f)();
    284354  }
    285355};
    286356template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)() const) const
    287357{
    288   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    289     ((*iter)->*f)();
     358  atom *Walker = start;
     359  while (Walker->next != end) {
     360    Walker = Walker->next;
     361    (Walker->*f)();
    290362  }
    291363};
     
    293365template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T), T t ) const
    294366{
    295   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    296     ((*iter)->*f)(t);
     367  atom *Walker = start;
     368  while (Walker->next != end) {
     369    Walker = Walker->next;
     370    (Walker->*f)(t);
    297371  }
    298372};
    299373template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T) const, T t ) const
    300374{
    301   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    302     ((*iter)->*f)(t);
     375  atom *Walker = start;
     376  while (Walker->next != end) {
     377    Walker = Walker->next;
     378    (Walker->*f)(t);
    303379  }
    304380};
     
    306382template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U), T t, U u ) const
    307383{
    308   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    309     ((*iter)->*f)(t, u);
     384  atom *Walker = start;
     385  while (Walker->next != end) {
     386    Walker = Walker->next;
     387    (Walker->*f)(t, u);
    310388  }
    311389};
    312390template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U) const, T t, U u ) const
    313391{
    314   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    315     ((*iter)->*f)(t, u);
     392  atom *Walker = start;
     393  while (Walker->next != end) {
     394    Walker = Walker->next;
     395    (Walker->*f)(t, u);
    316396  }
    317397};
     
    319399template <typename res, typename typ, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V), T t, U u, V v) const
    320400{
    321   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    322     ((*iter)->*f)(t, u, v);
     401  atom *Walker = start;
     402  while (Walker->next != end) {
     403    Walker = Walker->next;
     404    (Walker->*f)(t, u, v);
    323405  }
    324406};
    325407template <typename res, typename typ, typename T, typename U, typename V> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V) const, T t, U u, V v) const
    326408{
    327   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    328     ((*iter)->*f)(t, u, v);
     409  atom *Walker = start;
     410  while (Walker->next != end) {
     411    Walker = Walker->next;
     412    (Walker->*f)(t, u, v);
    329413  }
    330414};
     
    332416template <typename res, typename typ, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V, W), T t, U u, V v, W w) const
    333417{
    334   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    335     ((*iter)->*f)(t, u, v, w);
     418  atom *Walker = start;
     419  while (Walker->next != end) {
     420    Walker = Walker->next;
     421    (Walker->*f)(t, u, v, w);
    336422  }
    337423};
    338424template <typename res, typename typ, typename T, typename U, typename V, typename W> void molecule::ActOnAllAtoms( res (typ::*f)(T, U, V, W) const, T t, U u, V v, W w) const
    339425{
    340   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    341     ((*iter)->*f)(t, u, v, w);
     426  atom *Walker = start;
     427  while (Walker->next != end) {
     428    Walker = Walker->next;
     429    (Walker->*f)(t, u, v, w);
    342430  }
    343431};
     
    348436template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *) ) const
    349437{
     438  atom *Walker = start;
    350439  int inc = 1;
    351   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    352     (*Setor) (&array[((*iter)->*index)], &inc);
     440  while (Walker->next != end) {
     441    Walker = Walker->next;
     442    (*Setor) (&array[(Walker->*index)], &inc);
    353443  }
    354444};
    355445template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T value ) const
    356446{
    357   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    358     (*Setor) (&array[((*iter)->*index)], &value);
     447  atom *Walker = start;
     448  while (Walker->next != end) {
     449    Walker = Walker->next;
     450    (*Setor) (&array[(Walker->*index)], &value);
    359451  }
    360452};
    361453template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T *value ) const
    362454{
    363   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    364     (*Setor) (&array[((*iter)->*index)], value);
     455  atom *Walker = start;
     456  while (Walker->next != end) {
     457    Walker = Walker->next;
     458    (*Setor) (&array[(Walker->*index)], value);
    365459  }
    366460};
     
    368462template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *) ) const
    369463{
     464  atom *Walker = start;
    370465  int inc = 1;
    371   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    372     (*Setor) (&array[((*iter)->type->*index)], &inc);
     466  while (Walker->next != end) {
     467    Walker = Walker->next;
     468    (*Setor) (&array[(Walker->type->*index)], &inc);
    373469  }
    374470};
    375471template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T value ) const
    376472{
    377   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    378     (*Setor) (&array[((*iter)->type->*index)], &value);
     473  atom *Walker = start;
     474  while (Walker->next != end) {
     475    Walker = Walker->next;
     476    (*Setor) (&array[(Walker->type->*index)], &value);
    379477  }
    380478};
    381479template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T *value ) const
    382480{
    383   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    384     (*Setor) (&array[((*iter)->type->*index)], value);
     481  atom *Walker = start;
     482  while (Walker->next != end) {
     483    Walker = Walker->next;
     484    (*Setor) (&array[(Walker->type->*index)], value);
    385485  }
    386486};
     
    388488template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ atom::*value ) const
    389489{
    390   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    391     array[((*iter)->*index)] = ((*iter)->*Setor) ((*iter)->*value);
     490  atom *Walker = start;
     491  while (Walker->next != end) {
     492    Walker = Walker->next;
     493    array[(Walker->*index)] = (Walker->*Setor) (Walker->*value);
    392494  }
    393495};
    394496template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ atom::*value ) const
    395497{
    396   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    397     array[((*iter)->*index)] = ((*iter)->*Setor) ((*iter)->*value);
     498  atom *Walker = start;
     499  while (Walker->next != end) {
     500    Walker = Walker->next;
     501    array[(Walker->*index)] = (Walker->*Setor) (Walker->*value);
    398502  }
    399503};
    400504template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ &vect ) const
    401505{
    402   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    403     array[((*iter)->*index)] = ((*iter)->*Setor) (vect);
     506  atom *Walker = start;
     507  while (Walker->next != end) {
     508    Walker = Walker->next;
     509    array[(Walker->*index)] = (Walker->*Setor) (vect);
    404510  }
    405511};
    406512template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ &vect ) const
    407513{
    408   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    409     array[((*iter)->*index)] = ((*iter)->*Setor) (vect);
     514  atom *Walker = start;
     515  while (Walker->next != end) {
     516    Walker = Walker->next;
     517    array[(Walker->*index)] = (Walker->*Setor) (vect);
    410518  }
    411519};
    412520template <typename T, typename typ, typename typ2> void molecule::SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value ) const
    413521{
    414   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    415     (*iter)->*value = array[((*iter)->*index)];
    416     //Log() << Verbose(2) << *(*iter) << " gets " << ((*iter)->*value); << endl;
     522  atom *Walker = start;
     523  while (Walker->next != end) {
     524    Walker = Walker->next;
     525    Walker->*value = array[(Walker->*index)];
     526    //Log() << Verbose(2) << *Walker << " gets " << (Walker->*value); << endl;
    417527  }
    418528};
     
    420530template <typename T, typename typ> void molecule::SetAtomValueToValue ( T value, T typ::*ptr ) const
    421531{
    422   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    423     (*iter)->*ptr = value;
    424     //Log() << Verbose(2) << *(*iter) << " gets " << ((*iter)->*ptr) << endl;
     532  atom *Walker = start;
     533  while (Walker->next != end) {
     534    Walker = Walker->next;
     535    Walker->*ptr = value;
     536    //Log() << Verbose(2) << *Walker << " gets " << (Walker->*ptr) << endl;
    425537  }
    426538};
  • src/moleculelist.cpp

    ra7b761b r8f215d  
    2020#include "memoryallocator.hpp"
    2121#include "periodentafel.hpp"
    22 #include "Helpers/Assert.hpp"
     22#include "World.hpp"
    2323
    2424/*********************************** Functions for class MoleculeListClass *************************/
     
    6969  int Count, Counter, aCounter, bCounter;
    7070  int flag;
     71  atom *aWalker = NULL;
     72  atom *bWalker = NULL;
    7173
    7274  // sort each atom list and put the numbers into a list, then go through
    7375  //Log() << Verbose(0) << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
    74   // Yes those types are awkward... but check it for yourself it checks out this way
    75   molecule *const *mol1_ptr= static_cast<molecule *const *>(a);
    76   molecule *mol1 = *mol1_ptr;
    77   molecule *const *mol2_ptr= static_cast<molecule *const *>(b);
    78   molecule *mol2 = *mol2_ptr;
    79   if (mol1->getAtomCount() < mol2->getAtomCount()) {
     76  if ((**(molecule **) a).AtomCount < (**(molecule **) b).AtomCount) {
    8077    return -1;
    8178  } else {
    82     if (mol1->getAtomCount() > mol2->getAtomCount())
     79    if ((**(molecule **) a).AtomCount > (**(molecule **) b).AtomCount)
    8380      return +1;
    8481    else {
    85       Count = mol1->getAtomCount();
     82      Count = (**(molecule **) a).AtomCount;
    8683      aList = new int[Count];
    8784      bList = new int[Count];
    8885
    8986      // fill the lists
     87      aWalker = (**(molecule **) a).start;
     88      bWalker = (**(molecule **) b).start;
    9089      Counter = 0;
    9190      aCounter = 0;
    9291      bCounter = 0;
    93       molecule::const_iterator aiter = mol1->begin();
    94       molecule::const_iterator biter = mol2->begin();
    95       for (;(aiter != mol1->end()) && (biter != mol2->end());
    96           ++aiter, ++biter) {
    97         if ((*aiter)->GetTrueFather() == NULL)
     92      while ((aWalker->next != (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
     93        aWalker = aWalker->next;
     94        bWalker = bWalker->next;
     95        if (aWalker->GetTrueFather() == NULL)
    9896          aList[Counter] = Count + (aCounter++);
    9997        else
    100           aList[Counter] = (*aiter)->GetTrueFather()->nr;
    101         if ((*biter)->GetTrueFather() == NULL)
     98          aList[Counter] = aWalker->GetTrueFather()->nr;
     99        if (bWalker->GetTrueFather() == NULL)
    102100          bList[Counter] = Count + (bCounter++);
    103101        else
    104           bList[Counter] = (*biter)->GetTrueFather()->nr;
     102          bList[Counter] = bWalker->GetTrueFather()->nr;
    105103        Counter++;
    106104      }
    107105      // check if AtomCount was for real
    108106      flag = 0;
    109       if ((aiter == mol1->end()) && (biter != mol2->end())) {
     107      if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
    110108        flag = -1;
    111109      } else {
    112         if ((aiter != mol1->end()) && (biter == mol2->end()))
     110        if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **) b).end))
    113111          flag = 1;
    114112      }
     
    144142void MoleculeListClass::Enumerate(ostream *out)
    145143{
     144  atom *Walker = NULL;
    146145  periodentafel *periode = World::getInstance().getPeriode();
    147146  std::map<atomicNumber_t,unsigned int> counts;
     
    159158      // count atoms per element and determine size of bounding sphere
    160159      size=0.;
    161       for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
    162         counts[(*iter)->type->getNumber()]++;
    163         if ((*iter)->x.DistanceSquared(Origin) > size)
    164           size = (*iter)->x.DistanceSquared(Origin);
     160      Walker = (*ListRunner)->start;
     161      while (Walker->next != (*ListRunner)->end) {
     162        Walker = Walker->next;
     163        counts[Walker->type->getNumber()]++;
     164        if (Walker->x.DistanceSquared(Origin) > size)
     165          size = Walker->x.DistanceSquared(Origin);
    165166      }
    166167      // output Index, Name, number of atoms, chemical formula
    167       (*out) << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->getAtomCount() << "\t";
     168      (*out) << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->AtomCount << "\t";
    168169
    169170      std::map<atomicNumber_t,unsigned int>::reverse_iterator iter;
     
    201202
    202203  // put all molecules of src into mol
    203   molecule::iterator runner;
    204   for (molecule::iterator iter = srcmol->begin(); iter != srcmol->end(); ++iter) {
    205     runner = iter++;
    206     srcmol->UnlinkAtom((*runner));
    207     mol->AddAtom((*runner));
     204  atom *Walker = srcmol->start;
     205  atom *NextAtom = Walker->next;
     206  while (NextAtom != srcmol->end) {
     207    Walker = NextAtom;
     208    NextAtom = Walker->next;
     209    srcmol->UnlinkAtom(Walker);
     210    mol->AddAtom(Walker);
    208211  }
    209212
     
    225228
    226229  // put all molecules of src into mol
    227   atom *Walker = NULL;
    228   for (molecule::iterator iter = srcmol->begin(); iter != srcmol->end(); ++iter) {
    229     Walker = mol->AddCopyAtom((*iter));
     230  atom *Walker = srcmol->start;
     231  atom *NextAtom = Walker->next;
     232  while (NextAtom != srcmol->end) {
     233    Walker = NextAtom;
     234    NextAtom = Walker->next;
     235    Walker = mol->AddCopyAtom(Walker);
    230236    Walker->father = Walker;
    231237  }
     
    324330
    325331  // prepare index list for bonds
    326   atom ** CopyAtoms = new atom*[srcmol->getAtomCount()];
    327   for(int i=0;i<srcmol->getAtomCount();i++)
     332  srcmol->CountAtoms();
     333  atom ** CopyAtoms = new atom*[srcmol->AtomCount];
     334  for(int i=0;i<srcmol->AtomCount;i++)
    328335    CopyAtoms[i] = NULL;
    329336
    330337  // for each of the source atoms check whether we are in- or outside and add copy atom
     338  atom *Walker = srcmol->start;
    331339  int nr=0;
    332   for (molecule::const_iterator iter = srcmol->begin(); iter != srcmol->end(); ++iter) {
    333     DoLog(2) && (Log() << Verbose(2) << "INFO: Current Walker is " << **iter << "." << endl);
    334     if (!TesselStruct->IsInnerPoint((*iter)->x, LCList)) {
    335       CopyAtoms[(*iter)->nr] = (*iter)->clone();
    336       mol->AddAtom(CopyAtoms[(*iter)->nr]);
     340  while (Walker->next != srcmol->end) {
     341    Walker = Walker->next;
     342    DoLog(2) && (Log() << Verbose(2) << "INFO: Current Walker is " << *Walker << "." << endl);
     343    if (!TesselStruct->IsInnerPoint(Walker->x, LCList)) {
     344      CopyAtoms[Walker->nr] = Walker->clone();
     345      mol->AddAtom(CopyAtoms[Walker->nr]);
    337346      nr++;
    338347    } else {
     
    340349    }
    341350  }
    342   DoLog(1) && (Log() << Verbose(1) << nr << " of " << srcmol->getAtomCount() << " atoms have been merged.");
     351  DoLog(1) && (Log() << Verbose(1) << nr << " of " << srcmol->AtomCount << " atoms have been merged.");
    343352
    344353  // go through all bonds and add as well
     
    373382bool MoleculeListClass::AddHydrogenCorrection(char *path)
    374383{
     384  atom *Walker = NULL;
     385  atom *Runner = NULL;
    375386  bond *Binder = NULL;
    376387  double ***FitConstant = NULL, **correction = NULL;
     
    477488        correction[k][j] = 0.;
    478489    // 2. take every hydrogen that is a saturated one
    479     for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
    480       //Log() << Verbose(1) << "(*iter): " << *(*iter) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
    481       if (((*iter)->type->Z == 1) && (((*iter)->father == NULL)
    482           || ((*iter)->father->type->Z != 1))) { // if it's a hydrogen
    483         for (molecule::const_iterator runner = (*ListRunner)->begin(); runner != (*ListRunner)->end(); ++runner) {
    484           //Log() << Verbose(2) << "Runner: " << *(*runner) << " with first bond " << *((*iter)->ListOfBonds.begin()) << "." << endl;
     490    Walker = (*ListRunner)->start;
     491    while (Walker->next != (*ListRunner)->end) {
     492      Walker = Walker->next;
     493      //Log() << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
     494      if ((Walker->type->Z == 1) && ((Walker->father == NULL)
     495          || (Walker->father->type->Z != 1))) { // if it's a hydrogen
     496        Runner = (*ListRunner)->start;
     497        while (Runner->next != (*ListRunner)->end) {
     498          Runner = Runner->next;
     499          //Log() << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(Walker->ListOfBonds.begin()) << "." << endl;
    485500          // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
    486           Binder = *((*runner)->ListOfBonds.begin());
    487           if (((*runner)->type->Z == 1) && ((*runner)->nr > (*iter)->nr) && (Binder->GetOtherAtom((*runner)) != Binder->GetOtherAtom((*iter)))) { // (hydrogens have only one bonding partner!)
     501          Binder = *(Runner->ListOfBonds.begin());
     502          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (Binder->GetOtherAtom(Runner) != Binder->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
    488503            // 4. evaluate the morse potential for each matrix component and add up
    489             distance = (*runner)->x.distance((*iter)->x);
    490             //Log() << Verbose(0) << "Fragment " << (*ListRunner)->name << ": " << *(*runner) << "<= " << distance << "=>" << *(*iter) << ":" << endl;
     504            distance = Runner->x.distance(Walker->x);
     505            //Log() << Verbose(0) << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
    491506            for (int k = 0; k < a; k++) {
    492507              for (int j = 0; j < b; j++) {
     
    562577  ofstream ForcesFile;
    563578  stringstream line;
     579  atom *Walker = NULL;
    564580  periodentafel *periode=World::getInstance().getPeriode();
    565581
     
    574590      periodentafel::const_iterator elemIter;
    575591      for(elemIter=periode->begin();elemIter!=periode->end();++elemIter){
    576         if ((*ListRunner)->ElementsInMolecule[(*elemIter).first]) { // if this element got atoms
    577           for(molecule::iterator atomIter = (*ListRunner)->begin(); atomIter !=(*ListRunner)->end();++atomIter){
    578             if ((*atomIter)->type->getNumber() == (*elemIter).first) {
    579               if (((*atomIter)->GetTrueFather() != NULL) && ((*atomIter)->GetTrueFather() != (*atomIter))) {// if there is a rea
     592          if ((*ListRunner)->ElementsInMolecule[(*elemIter).first]) { // if this element got atoms
     593          Walker = (*ListRunner)->start;
     594          while (Walker->next != (*ListRunner)->end) { // go through every atom of this element
     595            Walker = Walker->next;
     596            if (Walker->type->getNumber() == (*elemIter).first) {
     597              if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
    580598                //Log() << Verbose(0) << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
    581                 ForcesFile << SortIndex[(*atomIter)->GetTrueFather()->nr] << "\t";
     599                ForcesFile << SortIndex[Walker->GetTrueFather()->nr] << "\t";
    582600              } else
    583601                // otherwise a -1 to indicate an added saturation hydrogen
     
    615633  bool result = true;
    616634  bool intermediateResult = true;
     635  atom *Walker = NULL;
    617636  Vector BoxDimension;
    618637  char *FragmentNumber = NULL;
     
    655674    // list atoms in fragment for debugging
    656675    DoLog(2) && (Log() << Verbose(2) << "Contained atoms: ");
    657     for (molecule::const_iterator iter = (*ListRunner)->begin(); iter != (*ListRunner)->end(); ++iter) {
    658       DoLog(0) && (Log() << Verbose(0) << (*iter)->getName() << " ");
     676    Walker = (*ListRunner)->start;
     677    while (Walker->next != (*ListRunner)->end) {
     678      Walker = Walker->next;
     679      DoLog(0) && (Log() << Verbose(0) << Walker->getName() << " ");
    659680    }
    660681    DoLog(0) && (Log() << Verbose(0) << endl);
     
    736757{
    737758  molecule *mol = World::getInstance().createMolecule();
     759  atom *Walker = NULL;
     760  atom *Advancer = NULL;
    738761  bond *Binder = NULL;
    739762  bond *Stepper = NULL;
     
    741764  for (MoleculeList::iterator MolRunner = ListOfMolecules.begin(); !ListOfMolecules.empty(); MolRunner = ListOfMolecules.begin()) {
    742765    // shift all atoms to new molecule
    743     for (molecule::iterator iter = (*MolRunner)->begin(); (*MolRunner)->empty(); iter = (*MolRunner)->begin()) {
    744       DoLog(3) && (Log() << Verbose(3) << "Re-linking " << (*iter) << "..." << endl);
    745       (*iter)->father = (*iter);
    746       mol->AddAtom((*iter));    // counting starts at 1
    747       (*MolRunner)->erase(iter);
     766    Advancer = (*MolRunner)->start->next;
     767    while (Advancer != (*MolRunner)->end) {
     768      Walker = Advancer;
     769      Advancer = Advancer->next;
     770      DoLog(3) && (Log() << Verbose(3) << "Re-linking " << *Walker << "..." << endl);
     771      unlink(Walker);
     772      Walker->father = Walker;
     773      mol->AddAtom(Walker);    // counting starts at 1
    748774    }
    749775    // remove all bonds
     
    799825  // 4b. create and fill map of which atom is associated to which connected molecule (note, counting starts at 1)
    800826  int FragmentCounter = 0;
    801   int *MolMap = Calloc<int>(mol->getAtomCount(), "config::Load() - *MolMap");
     827  int *MolMap = Calloc<int>(mol->AtomCount, "config::Load() - *MolMap");
    802828  MoleculeLeafClass *MolecularWalker = Subgraphs;
     829  Walker = NULL;
    803830  while (MolecularWalker->next != NULL) {
    804831    MolecularWalker = MolecularWalker->next;
    805     for (molecule::const_iterator iter = MolecularWalker->Leaf->begin(); iter != MolecularWalker->Leaf->end(); ++iter) {
    806       MolMap[(*iter)->GetTrueFather()->nr] = FragmentCounter+1;
     832    Walker = MolecularWalker->Leaf->start;
     833    while (Walker->next != MolecularWalker->Leaf->end) {
     834      Walker = Walker->next;
     835      MolMap[Walker->GetTrueFather()->nr] = FragmentCounter+1;
    807836    }
    808837    FragmentCounter++;
     
    810839
    811840  // 4c. relocate atoms to new molecules and remove from Leafs
    812   for (molecule::iterator iter = mol->begin(); !mol->empty(); iter = mol->begin()) {
    813     if (((*iter)->nr <0) || ((*iter)->nr >= mol->getAtomCount())) {
    814       DoeLog(0) && (eLog()<< Verbose(0) << "Index of atom " << **iter << " is invalid!" << endl);
     841  Walker = NULL;
     842  while (mol->start->next != mol->end) {
     843    Walker = mol->start->next;
     844    if ((Walker->nr <0) || (Walker->nr >= mol->AtomCount)) {
     845      DoeLog(0) && (eLog()<< Verbose(0) << "Index of atom " << *Walker << " is invalid!" << endl);
    815846      performCriticalExit();
    816847    }
    817     FragmentCounter = MolMap[(*iter)->nr];
     848    FragmentCounter = MolMap[Walker->nr];
    818849    if (FragmentCounter != 0) {
    819       DoLog(3) && (Log() << Verbose(3) << "Re-linking " << **iter << "..." << endl);
    820       molecules[FragmentCounter-1]->AddAtom((*iter));    // counting starts at 1
    821       mol->erase(iter);
     850      DoLog(3) && (Log() << Verbose(3) << "Re-linking " << *Walker << "..." << endl);
     851      unlink(Walker);
     852      molecules[FragmentCounter-1]->AddAtom(Walker);    // counting starts at 1
    822853    } else {
    823       DoeLog(0) && (eLog()<< Verbose(0) << "Atom " << **iter << " not associated to molecule!" << endl);
     854      DoeLog(0) && (eLog()<< Verbose(0) << "Atom " << *Walker << " not associated to molecule!" << endl);
    824855      performCriticalExit();
    825856    }
     
    829860  while (mol->first->next != mol->last) {
    830861    Binder = mol->first->next;
    831     const atom * const Walker = Binder->leftatom;
     862    Walker = Binder->leftatom;
    832863    unlink(Binder);
    833864    link(Binder,molecules[MolMap[Walker->nr]-1]->last);   // counting starts at 1
     
    851882int MoleculeListClass::CountAllAtoms() const
    852883{
     884  atom *Walker = NULL;
    853885  int AtomNo = 0;
    854886  for (MoleculeList::const_iterator MolWalker = ListOfMolecules.begin(); MolWalker != ListOfMolecules.end(); MolWalker++) {
    855     AtomNo += (*MolWalker)->size();
     887    Walker = (*MolWalker)->start;
     888    while (Walker->next != (*MolWalker)->end) {
     889      Walker = Walker->next;
     890      AtomNo++;
     891    }
    856892  }
    857893  return AtomNo;
     
    10281064bool MoleculeLeafClass::FillBondStructureFromReference(const molecule * const reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
    10291065{
     1066  atom *Walker = NULL;
    10301067  atom *OtherWalker = NULL;
    10311068  atom *Father = NULL;
     
    10351072  DoLog(1) && (Log() << Verbose(1) << "Begin of FillBondStructureFromReference." << endl);
    10361073  // fill ListOfLocalAtoms if NULL was given
    1037   if (!FillListOfLocalAtoms(ListOfLocalAtoms, FragmentCounter, reference->getAtomCount(), FreeList)) {
     1074  if (!FillListOfLocalAtoms(ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    10381075    DoLog(1) && (Log() << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl);
    10391076    return false;
     
    10511088    }
    10521089
    1053     for(molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end(); ++iter) {
    1054       Father = (*iter)->GetTrueFather();
     1090    Walker = Leaf->start;
     1091    while (Walker->next != Leaf->end) {
     1092      Walker = Walker->next;
     1093      Father = Walker->GetTrueFather();
    10551094      AtomNo = Father->nr; // global id of the current walker
    10561095      for (BondList::const_iterator Runner = Father->ListOfBonds.begin(); Runner != Father->ListOfBonds.end(); (++Runner)) {
    1057         OtherWalker = ListOfLocalAtoms[FragmentCounter][(*Runner)->GetOtherAtom((*iter)->GetTrueFather())->nr]; // local copy of current bond partner of walker
     1096        OtherWalker = ListOfLocalAtoms[FragmentCounter][(*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
    10581097        if (OtherWalker != NULL) {
    1059           if (OtherWalker->nr > (*iter)->nr)
    1060             Leaf->AddBond((*iter), OtherWalker, (*Runner)->BondDegree);
     1098          if (OtherWalker->nr > Walker->nr)
     1099            Leaf->AddBond(Walker, OtherWalker, (*Runner)->BondDegree);
    10611100        } else {
    1062           DoLog(1) && (Log() << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << (*Runner)->GetOtherAtom((*iter)->GetTrueFather())->nr << "] is NULL!" << endl);
     1101          DoLog(1) && (Log() << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << (*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl);
    10631102          status = false;
    10641103        }
     
    10871126bool MoleculeLeafClass::FillRootStackForSubgraphs(KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
    10881127{
    1089   atom *Father = NULL;
     1128  atom *Walker = NULL, *Father = NULL;
    10901129
    10911130  if (RootStack != NULL) {
     
    10931132    if (&(RootStack[FragmentCounter]) != NULL) {
    10941133      RootStack[FragmentCounter].clear();
    1095       for(molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end(); ++iter) {
    1096         Father = (*iter)->GetTrueFather();
     1134      Walker = Leaf->start;
     1135      while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
     1136        Walker = Walker->next;
     1137        Father = Walker->GetTrueFather();
    10971138        if (AtomMask[Father->nr]) // apply mask
    10981139#ifdef ADDHYDROGEN
    1099           if ((*iter)->type->Z != 1) // skip hydrogen
     1140          if (Walker->type->Z != 1) // skip hydrogen
    11001141#endif
    1101           RootStack[FragmentCounter].push_front((*iter)->nr);
     1142          RootStack[FragmentCounter].push_front(Walker->nr);
    11021143      }
    11031144      if (next != NULL)
     
    11381179
    11391180  if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
    1140     status = status && Leaf->CreateFatherLookupTable(ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
     1181    status = status && CreateFatherLookupTable(Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
    11411182    FreeList = FreeList && true;
    11421183  }
     
    11621203  DoLog(1) && (Log() << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl);
    11631204  // fill ListOfLocalAtoms if NULL was given
    1164   if (!FillListOfLocalAtoms(ListOfLocalAtoms, FragmentCounter, reference->getAtomCount(), FreeList)) {
     1205  if (!FillListOfLocalAtoms(ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
    11651206    DoLog(1) && (Log() << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl);
    11661207    return false;
  • src/tesselation.cpp

    ra7b761b r8f215d  
    2121#include "Plane.hpp"
    2222#include "Exceptions/LinearDependenceException.hpp"
    23 #include "Helpers/Assert.hpp"
    24 
    2523#include "Helpers/Assert.hpp"
    2624
     
    360358  // get ascending order of endpoints
    361359  PointMap OrderMap;
    362   for (int i = 0; i < 3; i++) {
     360  for (int i = 0; i < 3; i++)
    363361    // for all three lines
    364362    for (int j = 0; j < 2; j++) { // for both endpoints
     
    366364      // and we don't care whether insertion fails
    367365    }
    368   }
    369366  // set endpoints
    370367  int Counter = 0;
     
    375372    Counter++;
    376373  }
    377   ASSERT(Counter >= 3,"We have a triangle with only two distinct endpoints!");
    378 };
    379 
     374  if (Counter < 3) {
     375    DoeLog(0) && (eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl);
     376    performCriticalExit();
     377  }
     378}
     379;
    380380
    381381/** Destructor of BoundaryTriangleSet.
     
    12321232;
    12331233
     1234/** PointCloud implementation of GetTerminalPoint.
     1235 * Uses PointsOnBoundary and STL stuff.
     1236 */
     1237TesselPoint * Tesselation::GetTerminalPoint() const
     1238{
     1239  Info FunctionInfo(__func__);
     1240  PointMap::const_iterator Runner = PointsOnBoundary.end();
     1241  Runner--;
     1242  return (Runner->second->node);
     1243}
     1244;
     1245
    12341246/** PointCloud implementation of GoToNext.
    12351247 * Uses PointsOnBoundary and STL stuff.
     
    12431255;
    12441256
     1257/** PointCloud implementation of GoToPrevious.
     1258 * Uses PointsOnBoundary and STL stuff.
     1259 */
     1260void Tesselation::GoToPrevious() const
     1261{
     1262  Info FunctionInfo(__func__);
     1263  if (InternalPointer != PointsOnBoundary.begin())
     1264    InternalPointer--;
     1265}
     1266;
     1267
    12451268/** PointCloud implementation of GoToFirst.
    12461269 * Uses PointsOnBoundary and STL stuff.
     
    12501273  Info FunctionInfo(__func__);
    12511274  InternalPointer = PointsOnBoundary.begin();
     1275}
     1276;
     1277
     1278/** PointCloud implementation of GoToLast.
     1279 * Uses PointsOnBoundary and STL stuff.
     1280 */
     1281void Tesselation::GoToLast() const
     1282{
     1283  Info FunctionInfo(__func__);
     1284  InternalPointer = PointsOnBoundary.end();
     1285  InternalPointer--;
    12521286}
    12531287;
     
    25592593  // fill the set of neighbours
    25602594  TesselPointSet SetOfNeighbours;
    2561 
    25622595  SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
    25632596  for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
  • src/tesselation.hpp

    ra7b761b r8f215d  
    244244  virtual Vector *GetCenter() const { return NULL; };
    245245  virtual TesselPoint *GetPoint() const { return NULL; };
     246  virtual TesselPoint *GetTerminalPoint() const { return NULL; };
    246247  virtual int GetMaxId() const { return 0; };
    247248  virtual void GoToNext() const {};
     249  virtual void GoToPrevious() const {};
    248250  virtual void GoToFirst() const {};
     251  virtual void GoToLast() const {};
    249252  virtual bool IsEmpty() const { return true; };
    250253  virtual bool IsEnd() const { return true; };
     
    360363    virtual Vector *GetCenter(ofstream *out) const;
    361364    virtual TesselPoint *GetPoint() const;
     365    virtual TesselPoint *GetTerminalPoint() const;
    362366    virtual void GoToNext() const;
     367    virtual void GoToPrevious() const;
    363368    virtual void GoToFirst() const;
     369    virtual void GoToLast() const;
    364370    virtual bool IsEmpty() const;
    365371    virtual bool IsEnd() const;
  • src/tesselationhelpers.cpp

    ra7b761b r8f215d  
    838838    int i=cloud->GetMaxId();
    839839    int *LookupList = new int[i];
    840     for (cloud->GoToFirst(), i=0; !cloud->IsEnd(); cloud->GoToNext(), i++){
     840    for (cloud->GoToFirst(), i=0; !cloud->IsEnd(); cloud->GoToNext(), i++)
    841841      LookupList[i] = -1;
    842     }
    843842
    844843    // print atom coordinates
    845844    int Counter = 1;
    846845    TesselPoint *Walker = NULL;
    847     for (PointMap::const_iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); ++target) {
     846    for (PointMap::const_iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); target++) {
    848847      Walker = target->second->node;
    849848      LookupList[Walker->nr] = Counter++;
  • src/unittests/AnalysisCorrelationToPointUnitTest.cpp

    ra7b761b r8f215d  
    2525#include "periodentafel.hpp"
    2626#include "tesselation.hpp"
     27#include "World.hpp"
    2728
    2829#ifdef HAVE_TESTRUNNER
     
    7980
    8081  // check that TestMolecule was correctly constructed
    81   CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 4 );
     82  CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 );
    8283
    8384  TestList = World::getInstance().getMolecules();
  • src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

    ra7b761b r8f215d  
    2626#include "tesselation.hpp"
    2727#include "World.hpp"
    28 #include "Helpers/Assert.hpp"
    2928
    3029#include "Helpers/Assert.hpp"
     
    4140void AnalysisCorrelationToSurfaceUnitTest::setUp()
    4241{
    43   ASSERT_DO(Assert::Throw);
     42  //ASSERT_DO(Assert::Throw);
    4443
    4544  atom *Walker = NULL;
     
    7271  // construct molecule (tetraeder of hydrogens) base
    7372  TestSurfaceMolecule = World::getInstance().createMolecule();
    74 
    7573  Walker = World::getInstance().createAtom();
    7674  Walker->type = hydrogen;
    7775  *Walker->node = Vector(1., 0., 1. );
    78   TestSurfaceMolecule->AddAtom(Walker);
    79 
     76
     77  TestSurfaceMolecule->AddAtom(Walker);
    8078  Walker = World::getInstance().createAtom();
    8179  Walker->type = hydrogen;
    8280  *Walker->node = Vector(0., 1., 1. );
    8381  TestSurfaceMolecule->AddAtom(Walker);
    84 
    8582  Walker = World::getInstance().createAtom();
    8683  Walker->type = hydrogen;
    8784  *Walker->node = Vector(1., 1., 0. );
    8885  TestSurfaceMolecule->AddAtom(Walker);
    89 
    9086  Walker = World::getInstance().createAtom();
    9187  Walker->type = hydrogen;
     
    9490
    9591  // check that TestMolecule was correctly constructed
    96   CPPUNIT_ASSERT_EQUAL( TestSurfaceMolecule->getAtomCount(), 4 );
     92  CPPUNIT_ASSERT_EQUAL( TestSurfaceMolecule->AtomCount, 4 );
    9793
    9894  TestList = World::getInstance().getMolecules();
     
    111107  *Walker->node = Vector(4., 0., 4. );
    112108  TestSurfaceMolecule->AddAtom(Walker);
    113 
    114109  Walker = World::getInstance().createAtom();
    115110  Walker->type = carbon;
    116111  *Walker->node = Vector(0., 4., 4. );
    117112  TestSurfaceMolecule->AddAtom(Walker);
    118 
    119113  Walker = World::getInstance().createAtom();
    120114  Walker->type = carbon;
    121115  *Walker->node = Vector(4., 4., 0. );
    122116  TestSurfaceMolecule->AddAtom(Walker);
    123 
    124117  // add inner atoms
    125118  Walker = World::getInstance().createAtom();
     
    127120  *Walker->node = Vector(0.5, 0.5, 0.5 );
    128121  TestSurfaceMolecule->AddAtom(Walker);
    129 
    130122  TestSurfaceMolecule->ActiveFlag = true;
    131123  TestList->insert(TestSurfaceMolecule);
     
    158150void AnalysisCorrelationToSurfaceUnitTest::SurfaceTest()
    159151{
    160   CPPUNIT_ASSERT_EQUAL( 4, TestSurfaceMolecule->getAtomCount() );
     152  CPPUNIT_ASSERT_EQUAL( 4, TestSurfaceMolecule->AtomCount );
    161153  CPPUNIT_ASSERT_EQUAL( (size_t)2, TestList->ListOfMolecules.size() );
    162154  CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->PointsOnBoundary.size() );
  • src/unittests/AnalysisPairCorrelationUnitTest.cpp

    ra7b761b r8f215d  
    7979
    8080  // check that TestMolecule was correctly constructed
    81   CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 4 );
     81  CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 );
    8282
    8383  TestList = World::getInstance().getMolecules();
  • src/unittests/CountBondsUnitTest.cpp

    ra7b761b r8f215d  
    102102
    103103  // check that TestMolecule was correctly constructed
    104   CPPUNIT_ASSERT_EQUAL( TestMolecule1->getAtomCount(), 3 );
    105   CPPUNIT_ASSERT_EQUAL( TestMolecule2->getAtomCount(), 3 );
     104  CPPUNIT_ASSERT_EQUAL( TestMolecule1->AtomCount, 3 );
     105  Walker = TestMolecule1->start->next;
     106  CPPUNIT_ASSERT( TestMolecule1->end != Walker );
     107  CPPUNIT_ASSERT_EQUAL( TestMolecule2->AtomCount, 3 );
     108  Walker = TestMolecule2->start->next;
     109  CPPUNIT_ASSERT( TestMolecule2->end != Walker );
    106110
    107111  // create a small file with table
  • src/unittests/LinkedCellUnitTest.cpp

    ra7b761b r8f215d  
    6969
    7070  // check that TestMolecule was correctly constructed
    71   CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 3*3*3 );
     71  CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 3*3*3 );
     72  Walker = TestMolecule->start->next;
     73  CPPUNIT_ASSERT( TestMolecule->end != Walker );
    7274};
    7375
     
    195197{
    196198  // check all atoms
    197   for(molecule::iterator iter = TestMolecule->begin(); iter != TestMolecule->end();++iter){
    198     CPPUNIT_ASSERT_EQUAL( true, LC->SetIndexToNode(*iter) );
     199  atom *Walker = TestMolecule->start;
     200  while (Walker->next != TestMolecule->end) {
     201    Walker = Walker->next;
     202    CPPUNIT_ASSERT_EQUAL( true, LC->SetIndexToNode(Walker) );
    199203  }
    200204
    201205  // check internal vectors, returns false, because this atom is not in LC-list!
    202   atom *newAtom = World::getInstance().createAtom();
    203   newAtom->setName("test");
    204   newAtom->x= Vector(1,1,1);
    205   CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToNode(newAtom) );
    206   World::getInstance().destroyAtom(newAtom);
     206  Walker = World::getInstance().createAtom();
     207  Walker->setName("test");
     208  Walker->x= Vector(1,1,1);
     209  CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToNode(Walker) );
     210  World::getInstance().destroyAtom(Walker);
    207211
    208212  // check out of bounds vectors
    209   newAtom = World::getInstance().createAtom();
    210   newAtom->setName("test");
    211   newAtom->x = Vector(0,-1,0);
    212   CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToNode(newAtom) );
    213   World::getInstance().destroyAtom(newAtom);
     213  Walker = World::getInstance().createAtom();
     214  Walker->setName("test");
     215  Walker->x = Vector(0,-1,0);
     216  CPPUNIT_ASSERT_EQUAL( false, LC->SetIndexToNode(Walker) );
     217  World::getInstance().destroyAtom(Walker);
    214218};
    215219
     
    283287  size = ListOfPoints->size();
    284288  CPPUNIT_ASSERT_EQUAL( (size_t)27, size );
    285 
    286   for(molecule::iterator iter = TestMolecule->begin(); iter != TestMolecule->end(); ++iter){
    287     ListOfPoints->remove((*iter));
     289  Walker = TestMolecule->start;
     290  Walker = TestMolecule->start;
     291  while (Walker->next != TestMolecule->end) {
     292    Walker = Walker->next;
     293    ListOfPoints->remove(Walker);
    288294    size--;
    289295    CPPUNIT_ASSERT_EQUAL( size, ListOfPoints->size() );
     
    300306  size=ListOfPoints->size();
    301307  CPPUNIT_ASSERT_EQUAL( (size_t)8, size );
    302   for(molecule::iterator iter = TestMolecule->begin(); iter != TestMolecule->end(); ++iter){
    303     if (((*iter)->x[0] <2) && ((*iter)->x[1] <2) && ((*iter)->x[2] <2)) {
    304       ListOfPoints->remove(*iter);
     308  Walker = TestMolecule->start;
     309  while (Walker->next != TestMolecule->end) {
     310    Walker = Walker->next;
     311    if ((Walker->x[0] <2) && (Walker->x[1] <2) && (Walker->x[2] <2)) {
     312      ListOfPoints->remove(Walker);
    305313      size--;
    306314      CPPUNIT_ASSERT_EQUAL( size, ListOfPoints->size() );
     
    318326  size=ListOfPoints->size();
    319327  CPPUNIT_ASSERT_EQUAL( (size_t)27, size );
    320   for(molecule::iterator iter = TestMolecule->begin(); iter!=TestMolecule->end();++iter){
    321     ListOfPoints->remove(*iter);
     328  Walker = TestMolecule->start;
     329  while (Walker->next != TestMolecule->end) {
     330    Walker = Walker->next;
     331    ListOfPoints->remove(Walker);
    322332    size--;
    323333    CPPUNIT_ASSERT_EQUAL( size, ListOfPoints->size() );
     
    345355  size = ListOfPoints->size();
    346356  CPPUNIT_ASSERT_EQUAL( (size_t)7, size );
    347   for(molecule::iterator iter = TestMolecule->begin(); iter!=TestMolecule->end();++iter){
    348     if (((*iter)->x.DistanceSquared(tester) - 1.) < MYEPSILON ) {
    349       ListOfPoints->remove(*iter);
     357  Walker = TestMolecule->start;
     358  while (Walker->next != TestMolecule->end) {
     359    Walker = Walker->next;
     360    if ((Walker->x.DistanceSquared(tester) - 1.) < MYEPSILON ) {
     361      ListOfPoints->remove(Walker);
    350362      size--;
    351363      CPPUNIT_ASSERT_EQUAL( size, ListOfPoints->size() );
  • src/unittests/ObserverTest.cpp

    ra7b761b r8f215d  
    1111#include <cppunit/extensions/TestFactoryRegistry.h>
    1212#include <cppunit/ui/text/TestRunner.h>
    13 #include <set>
    1413
    1514#include "Patterns/Observer.hpp"
    16 #include "Patterns/ObservedIterator.hpp"
    1715#include "Helpers/Assert.hpp"
    1816
     
    164162  bool wasNotified;
    165163};
    166 
    167 class ObservableCollection : public Observable {
    168 public:
    169   typedef std::set<SimpleObservable*> set;
    170   typedef ObservedIterator<set> iterator;
    171   typedef set::const_iterator const_iterator;
    172 
    173   ObservableCollection(int _num) :
    174   num(_num)
    175   {
    176     for(int i=0; i<num; ++i){
    177       SimpleObservable *content = new SimpleObservable();
    178       content->signOn(this);
    179       theSet.insert(content);
    180     }
    181   }
    182 
    183   ~ObservableCollection(){
    184     set::iterator iter;
    185     for(iter=theSet.begin(); iter!=theSet.end(); ++iter ){
    186       delete (*iter);
    187     }
    188   }
    189 
    190   iterator begin(){
    191     return iterator(theSet.begin(),this);
    192   }
    193 
    194   iterator end(){
    195     return iterator(theSet.end(),this);
    196   }
    197 
    198   const int num;
    199 
    200 private:
    201   set theSet;
    202 };
    203 
    204164
    205165/******************* actuall tests ***************/
     
    213173  blockObservable = new BlockObservable();
    214174  notificationObservable = new NotificationObservable();
    215   collection = new ObservableCollection(5);
    216175
    217176  observer1 = new UpdateCountObserver();
     
    222181  notificationObserver1 = new NotificationObserver(notificationObservable->notification1);
    223182  notificationObserver2 = new NotificationObserver(notificationObservable->notification2);
     183
    224184}
    225185
     
    231191  delete blockObservable;
    232192  delete notificationObservable;
    233   delete collection;
    234193
    235194  delete observer1;
     
    318277  blockObservable->changeMethod2();
    319278  blockObservable->noChangeMethod();
    320 }
    321 
    322 void ObserverTest::iteratorTest(){
    323   int i = 0;
    324   // test the general iterator methods
    325   for(ObservableCollection::iterator iter=collection->begin(); iter!=collection->end();++iter){
    326     CPPUNIT_ASSERT(i< collection->num);
    327     i++;
    328   }
    329 
    330   i=0;
    331   for(ObservableCollection::const_iterator iter=collection->begin(); iter!=collection->end();++iter){
    332     CPPUNIT_ASSERT(i<collection->num);
    333     i++;
    334   }
    335 
    336   collection->signOn(observer1);
    337   {
    338     // we construct this out of the loop, so the iterator dies at the end of
    339     // the scope and not the end of the loop (allows more testing)
    340     ObservableCollection::iterator iter;
    341     for(iter=collection->begin(); iter!=collection->end(); ++iter){
    342       (*iter)->changeMethod();
    343     }
    344     // At this point no change should have been propagated
    345     CPPUNIT_ASSERT_EQUAL( 0, observer1->updates);
    346   }
    347   // After the Iterator has died the propagation should take place
    348   CPPUNIT_ASSERT_EQUAL( 1, observer1->updates);
    349 
    350   // when using a const_iterator no changes should be propagated
    351   for(ObservableCollection::const_iterator iter = collection->begin(); iter!=collection->end();++iter);
    352   CPPUNIT_ASSERT_EQUAL( 1, observer1->updates);
    353   collection->signOff(observer1);
    354279}
    355280
  • src/unittests/ObserverTest.hpp

    ra7b761b r8f215d  
    1717class CallObservable;
    1818class SuperObservable;
    19 class ObservableCollection;
    2019class BlockObservable;
    2120class NotificationObservable;
     21
    2222
    2323class ObserverTest :  public CppUnit::TestFixture
     
    2929  CPPUNIT_TEST ( doesNotifyTest );
    3030  CPPUNIT_TEST ( doesReportTest );
    31   CPPUNIT_TEST ( iteratorTest );
    3231  CPPUNIT_TEST ( CircleDetectionTest );
    3332  CPPUNIT_TEST_SUITE_END();
     
    4241  void doesNotifyTest();
    4342  void doesReportTest();
    44   void iteratorTest();
    4543  void CircleDetectionTest();
    4644
     
    6058  SuperObservable *superObservable;
    6159  NotificationObservable *notificationObservable;
    62   ObservableCollection *collection;
    63 
    6460};
    6561
  • src/unittests/analysisbondsunittest.cpp

    ra7b761b r8f215d  
    2525#include "molecule.hpp"
    2626#include "periodentafel.hpp"
    27 #include "World.hpp"
    2827
    2928#ifdef HAVE_TESTRUNNER
     
    9089
    9190  // check that TestMolecule was correctly constructed
    92   CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 5 );
     91  CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 5 );
    9392
    9493  // create a small file with table
  • src/unittests/bondgraphunittest.cpp

    ra7b761b r8f215d  
    8989
    9090  // check that TestMolecule was correctly constructed
    91   CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 4 );
     91  CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 );
    9292
    9393  // create a small file with table
     
    120120};
    121121
    122 /** Tests whether setup worked.
    123  */
    124 void BondGraphTest::SetupTest()
    125 {
    126   CPPUNIT_ASSERT_EQUAL (false, TestMolecule->empty());
    127   CPPUNIT_ASSERT_EQUAL ((size_t)4, TestMolecule->size());
    128 };
    129 
    130122/** UnitTest for BondGraphTest::LoadBondLengthTable().
    131123 */
     
    142134void BondGraphTest::ConstructGraphFromTableTest()
    143135{
    144   molecule::iterator Walker = TestMolecule->begin();
    145   molecule::iterator Runner = TestMolecule->begin();
    146   Runner++;
     136  atom *Walker = TestMolecule->start->next;
     137  atom *Runner = TestMolecule->end->previous;
     138  CPPUNIT_ASSERT( TestMolecule->end != Walker );
    147139  CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(*filename) );
    148140  CPPUNIT_ASSERT_EQUAL( true , BG->ConstructBondGraph(TestMolecule) );
    149   CPPUNIT_ASSERT_EQUAL( true , (*Walker)->IsBondedTo((*Runner)) );
     141  CPPUNIT_ASSERT_EQUAL( true , Walker->IsBondedTo(Runner) );
    150142};
    151143
     
    154146void BondGraphTest::ConstructGraphFromCovalentRadiiTest()
    155147{
    156 
    157   //atom *Walker = TestMolecule->start->next;
    158   //atom *Runner = TestMolecule->end->previous;
    159   //CPPUNIT_ASSERT( TestMolecule->end != Walker );
     148  atom *Walker = TestMolecule->start->next;
     149  atom *Runner = TestMolecule->end->previous;
     150  CPPUNIT_ASSERT( TestMolecule->end != Walker );
    160151  CPPUNIT_ASSERT_EQUAL( false , BG->LoadBondLengthTable(*dummyname) );
    161152  CPPUNIT_ASSERT_EQUAL( true , BG->ConstructBondGraph(TestMolecule) );
    162 
    163   // this cannot be assured using dynamic IDs
    164   //CPPUNIT_ASSERT_EQUAL( true , Walker->IsBondedTo(Runner) );
     153  CPPUNIT_ASSERT_EQUAL( true , Walker->IsBondedTo(Runner) );
    165154};
    166155
  • src/unittests/bondgraphunittest.hpp

    ra7b761b r8f215d  
    2222{
    2323    CPPUNIT_TEST_SUITE( BondGraphTest) ;
    24     CPPUNIT_TEST ( SetupTest );
    2524    CPPUNIT_TEST ( LoadTableTest );
    2625    CPPUNIT_TEST ( ConstructGraphFromTableTest );
     
    3130      void setUp();
    3231      void tearDown();
    33       void SetupTest();
    3432      void LoadTableTest();
    3533      void ConstructGraphFromTableTest();
  • src/unittests/listofbondsunittest.cpp

    ra7b761b r8f215d  
    7474
    7575  // check that TestMolecule was correctly constructed
    76   CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 4 );
     76  CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 );
    7777
    7878};
     
    9090};
    9191
    92 /** Tests whether setup worked correctly.
    93  *
    94  */
    95 void ListOfBondsTest::SetupTest()
    96 {
    97   CPPUNIT_ASSERT_EQUAL( false, TestMolecule->empty() );
    98   CPPUNIT_ASSERT_EQUAL( (size_t)4, TestMolecule->size() );
    99 };
    100 
    10192/** Unit Test of molecule::AddBond()
    10293 *
     
    10596{
    10697  bond *Binder = NULL;
    107   molecule::iterator iter = TestMolecule->begin();
    108   atom *atom1 = *iter;
    109   iter++;
    110   atom *atom2 = *iter;
     98  atom *atom1 = TestMolecule->start->next;
     99  atom *atom2 = atom1->next;
    111100  CPPUNIT_ASSERT( atom1 != NULL );
    112101  CPPUNIT_ASSERT( atom2 != NULL );
     
    135124{
    136125  bond *Binder = NULL;
    137   molecule::iterator iter = TestMolecule->begin();
    138   atom *atom1 = *iter;
    139   iter++;
    140   atom *atom2 = *iter;
     126  atom *atom1 = TestMolecule->start->next;
     127  atom *atom2 = atom1->next;
    141128  CPPUNIT_ASSERT( atom1 != NULL );
    142129  CPPUNIT_ASSERT( atom2 != NULL );
     
    163150{
    164151  bond *Binder = NULL;
    165   molecule::iterator iter = TestMolecule->begin();
    166   atom *atom1 = *iter;
    167   iter++;
    168   atom *atom2 = *iter;
    169   iter++;
    170   atom *atom3 = *iter;
     152  atom *atom1 = TestMolecule->start->next;
     153  atom *atom2 = atom1->next;
     154  atom *atom3 = atom2->next;
    171155  CPPUNIT_ASSERT( atom1 != NULL );
    172156  CPPUNIT_ASSERT( atom2 != NULL );
     
    205189{
    206190  bond *Binder = NULL;
    207   molecule::iterator iter = TestMolecule->begin();
    208   atom *atom1 = *iter;
    209   iter++;
    210   atom *atom2 = *iter;
     191  atom *atom1 = TestMolecule->start->next;
     192  atom *atom2 = atom1->next;
    211193  CPPUNIT_ASSERT( atom1 != NULL );
    212194  CPPUNIT_ASSERT( atom2 != NULL );
     
    233215{
    234216  bond *Binder = NULL;
    235   molecule::iterator iter = TestMolecule->begin();
    236   atom *atom1 = *iter;
    237   iter++;
    238   atom *atom2 = *iter;
     217  atom *atom1 = TestMolecule->start->next;
     218  atom *atom2 = atom1->next;
    239219  CPPUNIT_ASSERT( atom1 != NULL );
    240220  CPPUNIT_ASSERT( atom2 != NULL );
     
    260240{
    261241  bond *Binder = NULL;
    262   molecule::iterator iter = TestMolecule->begin();
    263   atom *atom1 = *iter;
    264   iter++;
    265   atom *atom2 = *iter;
     242  atom *atom1 = TestMolecule->start->next;
     243  atom *atom2 = atom1->next;
    266244  CPPUNIT_ASSERT( atom1 != NULL );
    267245  CPPUNIT_ASSERT( atom2 != NULL );
  • src/unittests/listofbondsunittest.hpp

    ra7b761b r8f215d  
    2020{
    2121    CPPUNIT_TEST_SUITE( ListOfBondsTest) ;
    22     CPPUNIT_TEST ( SetupTest );
    2322    CPPUNIT_TEST ( AddingBondTest );
    2423    CPPUNIT_TEST ( RemovingBondTest );
     
    3231      void setUp();
    3332      void tearDown();
    34       void SetupTest();
    3533      void AddingBondTest();
    3634      void RemovingBondTest();
Note: See TracChangeset for help on using the changeset viewer.