Changes in / [a7b761b:8f215d]
- Location:
- src
- Files:
-
- 1 added
- 1 deleted
- 39 edited
Legend:
- Unmodified
- Added
- Removed
-
src/Helpers/Assert.cpp
ra7b761b r8f215d 54 54 55 55 56 56 57 bool _my_assert::check(const bool res, 57 58 const char* condition, … … 62 63 { 63 64 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; 65 66 cout << "Assertion Message: " << message << std::endl; 66 67 while(true){ -
src/Legacy/oldmenu.cpp
ra7b761b r8f215d 423 423 void oldmenu::RemoveAtoms(molecule *mol) 424 424 { 425 atom * second;425 atom *first, *second; 426 426 int axis; 427 427 double tmp1, tmp2; … … 446 446 break; 447 447 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); 458 458 } 459 459 break; … … 465 465 Log() << Verbose(0) << "Upper boundary: "; 466 466 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); 473 475 } 474 476 } … … 513 515 min[i] = 0.; 514 516 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; 517 521 tmp1 = 0.; 518 if (first != (*iter)) {519 x = first->x - (*iter)->x;522 if (first != second) { 523 x = first->x - second->x; 520 524 tmp1 = x.Norm(); 521 525 } 522 526 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; 524 528 } 525 529 for (int i=MAX_ELEMENTS;i--;) … … 750 754 Log() << Verbose(0) << "State the factor: "; 751 755 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 754 760 Elements = new const element *[count]; 755 761 vectors = new Vector *[count]; 756 762 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; 760 768 j++; 761 769 } … … 1016 1024 return; 1017 1025 } 1026 atom *Walker = mol->start; 1018 1027 1019 1028 // generate some KeySets 1020 1029 Log() << Verbose(0) << "Generating KeySets." << endl; 1021 KeySet TestSets[mol-> getAtomCount()+1];1030 KeySet TestSets[mol->AtomCount+1]; 1022 1031 i=1; 1023 for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) { 1032 while (Walker->next != mol->end) { 1033 Walker = Walker->next; 1024 1034 for (int j=0;j<i;j++) { 1025 TestSets[j].insert( (*iter)->nr);1035 TestSets[j].insert(Walker->nr); 1026 1036 } 1027 1037 i++; … … 1029 1039 Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl; 1030 1040 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; 1039 1044 } 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); 1042 1049 1043 1050 // constructing Graph structure … … 1047 1054 // insert KeySets into Subgraphs 1048 1055 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++) { 1050 1057 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.))); 1051 1058 } 1052 1059 Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl; 1053 1060 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.))); 1055 1062 if (test2.second) { 1056 1063 Log() << Verbose(1) << "Insertion worked?!" << endl; -
src/Makefile.am
ra7b761b r8f215d 76 76 Descriptors/MoleculeIdDescriptor.cpp 77 77 78 78 79 79 DESCRIPTORHEADER = Descriptors/AtomDescriptor.hpp \ 80 80 Descriptors/AtomIdDescriptor.hpp \ … … 94 94 Exceptions/MathException.hpp \ 95 95 Exceptions/ZeroVectorException.hpp 96 97 96 98 97 SOURCE = ${ANALYSISSOURCE} \ … … 113 112 graph.cpp \ 114 113 helpers.cpp \ 115 Helpers/Assert.cpp \116 114 info.cpp \ 117 115 leastsquaremin.cpp \ 118 116 Line.cpp \ 119 117 linkedcell.cpp \ 118 lists.cpp \ 120 119 log.cpp \ 121 120 logger.cpp \ … … 134 133 tesselationhelpers.cpp \ 135 134 triangleintersectionlist.cpp \ 136 vector.cpp \137 135 verbose.cpp \ 138 136 vector_ops.cpp \ -
src/Patterns/Cacheable.hpp
ra7b761b r8f215d 28 28 owner(_owner) 29 29 {} 30 virtual T &getValue()=0;30 virtual T getValue()=0; 31 31 virtual void invalidate()=0; 32 32 virtual bool isValid()=0; … … 46 46 {} 47 47 48 virtual T &getValue(){48 virtual T getValue(){ 49 49 // set the state to valid 50 50 State::owner->switchState(State::owner->validState); … … 72 72 {} 73 73 74 virtual T &getValue(){74 virtual T getValue(){ 75 75 return content; 76 76 } … … 100 100 {} 101 101 102 virtual T &getValue(){102 virtual T getValue(){ 103 103 ASSERT(0,"Cannot get a value from a Cacheable after it's Observable has died"); 104 104 // we have to return a grossly invalid reference, because no value can be produced anymore … … 134 134 void subjectKilled(Observable *subject); 135 135 private: 136 136 137 void switchState(state_ptr newState); 137 138 … … 143 144 144 145 Observable *owner; 146 145 147 boost::function<T()> recalcMethod; 146 148 … … 219 221 220 222 const bool isValid() const; 221 const T &operator*() const;223 const T operator*() const; 222 224 223 225 // methods implemented for base-class Observer … … 235 237 236 238 template<typename T> 237 const T &Cacheable<T>::operator*() const{239 const T Cacheable<T>::operator*() const{ 238 240 return recalcMethod(); 239 241 } -
src/Patterns/Observer.cpp
ra7b761b r8f215d 82 82 Observable::_Observable_protector::_Observable_protector(Observable *_protege) : 83 83 protege(_protege) 84 {85 start_observer_internal(protege);86 }87 88 Observable::_Observable_protector::_Observable_protector(const _Observable_protector &dest) :89 protege(dest.protege)90 84 { 91 85 start_observer_internal(protege); … … 195 189 ASSERT(callTable.count(this),"SignOff called for an Observable without Observers."); 196 190 callees_t &callees = callTable[this]; 197 198 191 callees_t::iterator iter; 199 192 callees_t::iterator deliter; -
src/Patterns/Observer.hpp
ra7b761b r8f215d 36 36 typedef Notification *const Notification_ptr; 37 37 38 template<class _Set>39 class ObservedIterator;40 41 38 /** 42 39 * An Observer is notified by all Observed objects, when anything changes. … … 56 53 friend class Observable; 57 54 friend class Notification; 58 template<class> friend class ObservedIterator;59 60 55 public: 61 56 Observer(); … … 157 152 static std::set<Observable*> busyObservables; 158 153 154 159 155 //! @cond 160 156 // Structure for RAII-Style notification … … 168 164 public: 169 165 _Observable_protector(Observable *); 170 _Observable_protector(const _Observable_protector&);171 166 ~_Observable_protector(); 172 167 private: -
src/analysis_bonds.cpp
ra7b761b r8f215d 26 26 Mean = 0.; 27 27 28 atom *Walker = mol->start; 28 29 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(); 31 33 if (Max < count) 32 34 Max = count; … … 56 58 57 59 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) { 62 66 const double distance = (*BondRunner)->GetDistanceSquared(); 63 67 if (Min > distance) … … 122 126 int CountHydrogenBridgeBonds(MoleculeListClass *molecules, element * InterfaceElement = NULL) 123 127 { 128 atom *Walker = NULL; 129 atom *Runner = NULL; 124 130 int count = 0; 125 131 int OtherHydrogens = 0; … … 127 133 bool InterfaceFlag = false; 128 134 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)) { 136 144 // check distance 137 const double distance = (*Runner)->x.DistanceSquared((*Walker)->x);145 const double distance = Runner->x.DistanceSquared(Walker->x); 138 146 if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means different atoms 139 147 // on other atom(Runner) we check for bond to interface element and … … 143 151 OtherHydrogens = 0; 144 152 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); 147 155 // if hydrogen, check angle to be greater(!) than 30 degrees 148 156 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); 150 158 OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON); 151 159 Otherangle += angle; … … 168 176 if (InterfaceFlag && OtherHydrogenFlag) { 169 177 // 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); 172 180 if (OtherAtom->type->Z == 1) { 173 181 // 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); 176 184 count++; 177 185 break; … … 197 205 int CountBondsOfTwo(MoleculeListClass * const molecules, const element * const first, const element * const second) 198 206 { 207 atom *Walker = NULL; 199 208 int count = 0; 200 209 201 210 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 matches206 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)) { 209 218 count++; 210 219 DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << " bond found between " << *Walker << " and " << *OtherAtom << "." << endl); … … 231 240 bool MatchFlag[2]; 232 241 bool result = false; 242 atom *Walker = NULL; 233 243 const element * ElementArray[2]; 234 244 ElementArray[0] = first; … … 236 246 237 247 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 matches248 Walker = (*MolWalker)->start; 249 while (Walker->next != (*MolWalker)->end) { 250 Walker = Walker->next; 251 if (Walker->type == second) { // first element matches 242 252 for (int i=0;i<2;i++) 243 253 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); 246 256 for (int i=0;i<2;i++) 247 257 if ((!MatchFlag[i]) && (OtherAtom->type == ElementArray[i])) { -
src/analysis_correlation.cpp
ra7b761b r8f215d 40 40 } 41 41 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++) 43 43 if ((*MolWalker)->ActiveFlag) { 44 44 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++) 50 51 if ((*MolOtherWalker)->ActiveFlag) { 51 52 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) ) ); 59 62 } 60 }61 63 } 62 }63 64 } 64 65 } 65 66 } 66 67 } 67 } 68 68 69 return outmap; 69 70 }; … … 100 101 double * FullInverseMatrix = InverseMatrix(FullMatrix); 101 102 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); 106 109 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3 107 110 // go through every range in xyz and get distance … … 114 117 if ((*MolOtherWalker)->ActiveFlag) { 115 118 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); 121 126 periodicOtherX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3 122 127 // go through every range in xyz and get distance … … 127 132 checkOtherX.MatrixMultiplication(FullMatrix); 128 133 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) ) ); 131 136 } 132 137 } … … 164 169 if ((*MolWalker)->ActiveFlag) { 165 170 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()); 170 177 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) ) ); 172 179 } 173 180 } … … 204 211 double * FullInverseMatrix = InverseMatrix(FullMatrix); 205 212 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); 210 219 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3 211 220 // go through every range in xyz and get distance … … 217 226 distance = checkX.distance(*point); 218 227 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) ) ); 220 229 } 221 230 } … … 252 261 if ((*MolWalker)->ActiveFlag) { 253 262 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); 258 269 distance = Intersections.GetSmallestDistance(); 259 270 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) ) ); 261 272 } 262 273 } … … 303 314 double * FullInverseMatrix = InverseMatrix(FullMatrix); 304 315 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); 309 322 periodicX.MatrixMultiplication(FullInverseMatrix); // x now in [0,1)^3 310 323 // go through every range in xyz and get distance … … 324 337 } 325 338 // 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) ) ); 327 340 //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl; 328 341 } -
src/atom.cpp
ra7b761b r8f215d 59 59 atom::~atom() 60 60 { 61 unlink(this); 61 62 }; 62 63 -
src/bondgraph.cpp
ra7b761b r8f215d 90 90 bool status = true; 91 91 92 if (mol-> empty()) // only construct if molecule is not empty92 if (mol->start->next == mol->end) // only construct if molecule is not empty 93 93 return false; 94 94 … … 124 124 max_distance = 0.; 125 125 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; 129 131 } 130 132 max_distance *= 2.; -
src/boundary.cpp
ra7b761b r8f215d 139 139 { 140 140 Info FunctionInfo(__func__); 141 atom *Walker = NULL; 141 142 PointMap PointsOnBoundary; 142 143 LineMap LinesOnBoundary; … … 164 165 165 166 // 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); 168 171 ProjectedVector.ProjectOntoPlane(AxisVector); 169 172 … … 179 182 angle = 2. * M_PI - angle; 180 183 } 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))); 183 186 if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity 184 187 DoLog(2) && (Log() << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl); 185 188 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); 187 190 const double ProjectedVectorNorm = ProjectedVector.NormSquared(); 188 191 if ((ProjectedVectorNorm - BoundaryTestPair.first->second.first) > MYEPSILON) { 189 192 BoundaryTestPair.first->second.first = ProjectedVectorNorm; 190 BoundaryTestPair.first->second.second = (*iter);193 BoundaryTestPair.first->second.second = Walker; 191 194 DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << "." << endl); 192 195 } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) { 193 helper = (*iter)->x; 194 helper -= *MolCenter; 196 helper = Walker->x - (*MolCenter); 195 197 const double oldhelperNorm = helper.NormSquared(); 196 198 helper = BoundaryTestPair.first->second.second->x - (*MolCenter); 197 199 if (helper.NormSquared() < oldhelperNorm) { 198 BoundaryTestPair.first->second.second = (*iter);200 BoundaryTestPair.first->second.second = Walker; 199 201 DoLog(2) && (Log() << Verbose(2) << "Keeping new vector due to larger distance to molecule center " << helper.NormSquared() << "." << endl); 200 202 } else { … … 693 695 int repetition[NDIM] = { 1, 1, 1 }; 694 696 int TotalNoClusters = 1; 697 atom *Walker = NULL; 695 698 double totalmass = 0.; 696 699 double clustervolume = 0.; … … 716 719 717 720 // 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; 720 725 } 721 726 DoLog(0) && (Log() << Verbose(0) << "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit." << endl); … … 799 804 Vector Inserter; 800 805 double FillIt = false; 806 atom *Walker = NULL; 801 807 bond *Binder = NULL; 802 808 double phi[NDIM]; … … 805 811 806 812 for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) 807 if ((*ListRunner)-> getAtomCount()> 0) {813 if ((*ListRunner)->AtomCount > 0) { 808 814 DoLog(1) && (Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl); 809 815 LCList[(*ListRunner)] = new LinkedCell((*ListRunner), 10.); // get linked cell list … … 823 829 } 824 830 825 atom * CopyAtoms[filler->getAtomCount()]; 831 filler->CountAtoms(); 832 atom * CopyAtoms[filler->AtomCount]; 826 833 827 834 // calculate filler grid in [0,1]^3 … … 848 855 849 856 // go through all atoms 850 for (int i=0;i<filler-> getAtomCount();i++)857 for (int i=0;i<filler->AtomCount;i++) 851 858 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; 853 862 854 863 // create atomic random translation vector ... … … 874 883 875 884 // ... and put at new position 876 Inserter = (*iter)->x;885 Inserter = Walker->x; 877 886 if (DoRandomRotation) 878 887 Inserter.MatrixMultiplication(Rotations); … … 898 907 DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is outer point." << endl); 899 908 // 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); 904 913 } else { 905 914 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; 907 916 continue; 908 917 } … … 943 952 bool TesselationFailFlag = false; 944 953 945 mol->getAtomCount();946 947 954 if (TesselStruct == NULL) { 948 955 DoLog(1) && (Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl); … … 1016 1023 // // look whether all points are inside of the convex envelope, otherwise add them via degenerated triangles 1017 1024 // //->InsertStraddlingPoints(mol, LCList); 1018 // for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {1025 // mol->GoToFirst(); 1019 1026 // class TesselPoint *Runner = NULL; 1020 // Runner = *iter; 1027 // while (!mol->IsEnd()) { 1028 // Runner = mol->GetPoint(); 1021 1029 // Log() << Verbose(1) << "Checking on " << Runner->Name << " ... " << endl; 1022 1030 // if (!->IsInnerPoint(Runner, LCList)) { … … 1026 1034 // Log() << Verbose(2) << Runner->Name << " is inside of or on envelope." << endl; 1027 1035 // } 1036 // mol->GoToNext(); 1028 1037 // } 1029 1038 … … 1034 1043 status = CheckListOfBaselines(TesselStruct); 1035 1044 1036 cout << "before correction" << endl;1037 1038 1045 // store before correction 1039 StoreTrianglesinFile(mol, TesselStruct, filename, "");1046 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, ""); 1040 1047 1041 1048 // // correct degenerated polygons … … 1047 1054 // write final envelope 1048 1055 CalculateConcavityPerBoundaryPoint(TesselStruct); 1049 cout << "after correction" << endl; 1050 StoreTrianglesinFile(mol, TesselStruct, filename, ""); 1056 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, ""); 1051 1057 1052 1058 if (freeLC) -
src/builder.cpp
ra7b761b r8f215d 1744 1744 if (first->type != NULL) { 1745 1745 mol->AddAtom(first); // add to molecule 1746 if ((configPresent == empty) && (mol-> getAtomCount()!= 0))1746 if ((configPresent == empty) && (mol->AtomCount != 0)) 1747 1747 configPresent = present; 1748 1748 } else … … 1773 1773 DoLog(1) && (Log() << Verbose(1) << "Depth-First-Search Analysis." << endl); 1774 1774 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis 1775 int *MinimumRingSize = new int[mol-> getAtomCount()];1775 int *MinimumRingSize = new int[mol->AtomCount]; 1776 1776 atom ***ListOfLocalAtoms = NULL; 1777 1777 class StackClass<bond *> *BackEdgeStack = NULL; … … 1921 1921 int counter = 0; 1922 1922 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)) { 1924 1924 Boundary = *BigFinder; 1925 1925 } … … 1980 1980 performCriticalExit(); 1981 1981 } else { 1982 mol->getAtomCount();1983 1982 SaveFlag = true; 1984 1983 DoLog(1) && (Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl); … … 2099 2098 int counter = 0; 2100 2099 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)) { 2102 2102 Boundary = *BigFinder; 2103 2103 } 2104 2104 counter++; 2105 2105 } 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); 2107 2107 start = clock(); 2108 2108 LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.); … … 2180 2180 double tmp1 = atof(argv[argptr+1]); 2181 2181 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); 2190 2190 } 2191 2191 } else { … … 2332 2332 performCriticalExit(); 2333 2333 } else { 2334 mol->getAtomCount();2335 2334 SaveFlag = true; 2336 2335 DoLog(1) && (Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl); … … 2452 2451 faktor = 1; 2453 2452 } 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 2456 2456 Elements = new const element *[count]; 2457 2457 vectors = new Vector *[count]; 2458 2458 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; 2462 2464 j++; 2463 2465 } … … 2554 2556 int main(int argc, char **argv) 2555 2557 { 2556 // while we are non interactive, we want to abort from asserts2557 //ASSERT_DO(Assert::Abort);2558 2558 molecule *mol = NULL; 2559 2559 config *configuration = new config; … … 2595 2595 } 2596 2596 2597 // In the interactive mode, we can leave the user the choice in case of error2598 ASSERT_DO(Assert::Ask);2599 2600 2597 { 2601 2598 cout << ESPACKVersion << endl; -
src/config.cpp
ra7b761b r8f215d 1547 1547 int AtomNo = -1; 1548 1548 int MolNo = 0; 1549 atom *Walker = NULL; 1549 1550 FILE *f = NULL; 1550 1551 … … 1559 1560 fprintf(f, "# Created by MoleCuilder\n"); 1560 1561 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; 1562 1564 int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo"); 1563 1565 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 1567 1570 fprintf(f, 1568 1571 "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 */ 1570 1573 name, /* atom name */ 1571 (* MolRunner)->name, /* residue name */1574 (*Runner)->name, /* residue name */ 1572 1575 'a'+(unsigned char)(AtomNo % 26), /* letter for chain */ 1573 1576 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 */ 1579 1582 "0", /* segment identifier */ 1580 (*iter)->type->symbol, /* element symbol */1583 Walker->type->symbol, /* element symbol */ 1581 1584 "0"); /* charge */ 1582 1585 AtomNo++; … … 1598 1601 { 1599 1602 int AtomNo = -1; 1603 atom *Walker = NULL; 1600 1604 FILE *f = NULL; 1601 1605 … … 1612 1616 fprintf(f, "# Created by MoleCuilder\n"); 1613 1617 1618 Walker = mol->start; 1614 1619 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 1618 1624 fprintf(f, 1619 1625 "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 */ 1621 1627 name, /* atom name */ 1622 1628 mol->name, /* residue name */ 1623 1629 'a'+(unsigned char)(AtomNo % 26), /* letter for chain */ 1624 1630 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 */ 1630 1636 "0", /* segment identifier */ 1631 (*iter)->type->symbol, /* element symbol */1637 Walker->type->symbol, /* element symbol */ 1632 1638 "0"); /* charge */ 1633 1639 AtomNo++; … … 1647 1653 bool config::SaveTREMOLO(const char * const filename, const molecule * const mol) const 1648 1654 { 1655 atom *Walker = NULL; 1649 1656 ofstream *output = NULL; 1650 1657 stringstream * const fname = new stringstream; … … 1659 1666 1660 1667 // scan maximum number of neighbours 1668 Walker = mol->start; 1661 1669 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(); 1664 1673 if (MaxNeighbours < count) 1665 1674 MaxNeighbours = count; 1666 1675 } 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"; 1672 1683 *output << mol->name << "\t"; 1673 1684 *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++) 1680 1691 *output << "-\t"; 1681 1692 *output << endl; … … 1697 1708 bool config::SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const 1698 1709 { 1710 atom *Walker = NULL; 1699 1711 ofstream *output = NULL; 1700 1712 stringstream * const fname = new stringstream; … … 1711 1723 int MaxNeighbours = 0; 1712 1724 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(); 1715 1729 if (MaxNeighbours < count) 1716 1730 MaxNeighbours = count; 1717 1731 } 1718 1732 } 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; 1720 1734 1721 1735 // create global to local id map … … 1738 1752 int AtomNo = 0; 1739 1753 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; 1741 1757 *output << AtomNo+1 << "\t"; 1742 *output << (*iter)->getName() << "\t";1758 *output << Walker->getName() << "\t"; 1743 1759 *output << (*MolWalker)->name << "\t"; 1744 1760 *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++) 1751 1767 *output << "-\t"; 1752 1768 *output << endl; -
src/helpers.hpp
ra7b761b r8f215d 83 83 }; 84 84 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 */ 92 template <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 85 141 /** Frees a two-dimensional array. 86 142 * \param *ptr pointer to array -
src/lists.hpp
ra7b761b r8f215d 134 134 }; 135 135 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 */ 140 template <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 */ 152 template <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 136 160 #endif /* LISTS_HPP_ */ -
src/molecule.cpp
ra7b761b r8f215d 35 35 * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero. 36 36 */ 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),37 molecule::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), 39 39 BondCount(0), ElementCount(0), NoNonHydrogen(0), NoNonBonds(0), NoCyclicBonds(0), BondDistance(0.), 40 40 ActiveFlag(false), IndexNr(-1), 41 41 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 44 50 // init bond chain list 45 51 link(first,last); … … 63 69 delete(first); 64 70 delete(last); 71 end->getWorld()->destroyAtom(end); 72 start->getWorld()->destroyAtom(start); 65 73 }; 66 74 … … 75 83 } 76 84 77 int molecule::getAtomCount() const{78 return *AtomCount;79 }80 81 85 void molecule::setName(const std::string _name){ 82 86 OBSERVE; … … 100 104 stringstream sstr; 101 105 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()]++; 104 108 } 105 109 std::map<atomicNumber_t,unsigned int>::reverse_iterator iter; … … 111 115 } 112 116 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() const133 {134 return (begin() == end());135 }136 137 size_t molecule::size() const138 {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 ) const165 {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 }174 117 175 118 /** Adds given atom \a *pointer from molecule list. … … 180 123 bool molecule::AddAtom(atom *pointer) 181 124 { 125 bool retval = false; 182 126 OBSERVE; 183 127 if (pointer != NULL) { 184 128 pointer->sort = &pointer->nr; 129 pointer->nr = last_atom++; // increase number within molecule 130 AtomCount++; 185 131 if (pointer->type != NULL) { 186 132 if (ElementsInMolecule[pointer->type->Z] == 0) … … 195 141 } 196 142 } 197 insert(pointer);198 } 199 return true;143 retval = add(pointer, end); 144 } 145 return retval; 200 146 }; 201 147 … … 211 157 if (pointer != NULL) { 212 158 atom *walker = pointer->clone(); 213 walker->setName(pointer->getName()); 159 stringstream sstr; 160 sstr << pointer->getName(); 161 walker->setName(sstr.str()); 214 162 walker->nr = last_atom++; // increase number within molecule 215 insert(walker);163 add(walker, end); 216 164 if ((pointer->type != NULL) && (pointer->type->Z != 1)) 217 165 NoNonHydrogen++; 166 AtomCount++; 218 167 retval=walker; 219 168 } … … 625 574 626 575 // copy values 576 copy->CountAtoms(); 627 577 copy->CountElements(); 628 578 if (first->next != last) { // if adjaceny list is present … … 659 609 { 660 610 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 } 675 621 return Binder; 676 622 }; … … 746 692 bool molecule::RemoveAtom(atom *pointer) 747 693 { 748 ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");749 OBSERVE;750 694 if (ElementsInMolecule[pointer->type->Z] != 0) { // this would indicate an error 751 695 ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element 696 AtomCount--; 752 697 } else 753 698 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); … … 755 700 ElementCount--; 756 701 RemoveBonds(pointer); 757 erase(pointer); 758 return true; 702 return remove(pointer, start, end); 759 703 }; 760 704 … … 773 717 if (ElementsInMolecule[pointer->type->Z] == 0) // was last atom of this element? 774 718 ElementCount--; 775 erase(pointer);719 unlink(pointer); 776 720 return true; 777 721 }; … … 782 726 bool molecule::CleanupMolecule() 783 727 { 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)); 787 729 }; 788 730 … … 791 733 * \return pointer to atom or NULL 792 734 */ 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()) { 735 atom * molecule::FindAtom(int Nr) const{ 736 atom * walker = find(&Nr, start,end); 737 if (walker != NULL) { 800 738 //Log() << Verbose(0) << "Found Atom Nr. " << walker->nr << endl; 801 return (*iter);739 return walker; 802 740 } else { 803 741 DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl); … … 929 867 now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time' 930 868 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); 932 870 ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step ); 933 871 } … … 946 884 if (output != NULL) { 947 885 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); 949 887 ActOnAllAtoms( &atom::OutputXYZLine, output ); 950 888 return true; … … 956 894 * \param *out output stream for debugging 957 895 */ 958 int molecule::doCountAtoms() 959 { 960 int res = size(); 896 void molecule::CountAtoms() 897 { 961 898 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; 971 902 i++; 972 903 } 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 } 974 927 }; 975 928 … … 1033 986 /// first count both their atoms and elements and update lists thereby ... 1034 987 //Log() << Verbose(0) << "Counting atoms, updating list" << endl; 988 CountAtoms(); 989 OtherMolecule->CountAtoms(); 1035 990 CountElements(); 1036 991 OtherMolecule->CountElements(); … … 1039 994 /// -# AtomCount 1040 995 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); 1043 998 result = false; 1044 } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount()<< endl;999 } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl; 1045 1000 } 1046 1001 /// -# ElementCount … … 1079 1034 if (result) { 1080 1035 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"); 1083 1038 SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity); 1084 1039 SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity); … … 1086 1041 /// ... sort each list (using heapsort (o(N log N)) from GSL) 1087 1042 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"); 1093 1048 DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl); 1094 for(int i= getAtomCount();i--;)1049 for(int i=AtomCount;i--;) 1095 1050 PermutationMap[PermMap[i]] = (int) OtherPermMap[i]; 1096 1051 … … 1098 1053 DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl); 1099 1054 flag = 0; 1100 for (int i=0;i< getAtomCount();i++) {1055 for (int i=0;i<AtomCount;i++) { 1101 1056 DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl); 1102 1057 if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold) … … 1134 1089 int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule) 1135 1090 { 1091 atom *Walker = NULL, *OtherWalker = NULL; 1136 1092 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--;) 1139 1095 AtomicMap[i] = -1; 1140 1096 if (OtherMolecule == this) { // same molecule 1141 for (int i= getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence1097 for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence 1142 1098 AtomicMap[i] = i; 1143 1099 DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl); 1144 1100 } else { 1145 1101 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; 1149 1107 } 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; 1151 1111 //for (int i=0;i<AtomCount;i++) { // search atom 1152 1112 //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; 1156 1116 } 1157 1117 } 1158 DoLog(0) && (Log() << Verbose(0) << AtomicMap[ (*iter)->nr] << "\t");1118 DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t"); 1159 1119 } 1160 1120 DoLog(0) && (Log() << Verbose(0) << endl); … … 1190 1150 void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const 1191 1151 { 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; 1194 1156 } 1195 1157 }; -
src/molecule.hpp
ra7b761b r8f215d 34 34 #include "tesselation.hpp" 35 35 #include "Patterns/Observer.hpp" 36 #include "Patterns/ObservedIterator.hpp"37 36 #include "Patterns/Cacheable.hpp" 38 37 … … 89 88 friend molecule *NewMolecule(); 90 89 friend void DeleteMolecule(molecule *); 91 92 90 public: 93 typedef std::set<atom*> atomSet;94 typedef ObservedIterator<atomSet> iterator;95 typedef atomSet::const_iterator const_iterator;96 97 91 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 101 94 bond *first; //!< start of bond list 102 95 bond *last; //!< end of bond list 103 96 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() 105 98 int BondCount; //!< number of atoms, brought up-to-date by CountBonds() 106 99 int ElementCount; //!< how many unique elements are therein … … 117 110 private: 118 111 Cacheable<string> formula; 119 Cacheable<int> AtomCount;120 112 moleculeId_t id; 121 atomSet atoms; //<!set of atoms122 113 protected: 123 //void CountAtoms();124 /**125 * this iterator type should be used for internal variables, \126 * since it will not lock127 */128 typedef atomSet::iterator internal_iterator;129 130 131 114 molecule(const periodentafel * const teil); 132 115 virtual ~molecule(); … … 136 119 //getter and setter 137 120 const std::string getName(); 138 int getAtomCount() const;139 int doCountAtoms();140 121 moleculeId_t getId(); 141 122 void setId(moleculeId_t); … … 144 125 std::string calcFormula(); 145 126 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 157 127 158 128 // re-definition of virtual functions from PointCloud … … 160 130 Vector *GetCenter() const ; 161 131 TesselPoint *GetPoint() const ; 132 TesselPoint *GetTerminalPoint() const ; 162 133 int GetMaxId() const; 163 134 void GoToNext() const ; 135 void GoToPrevious() const ; 164 136 void GoToFirst() const ; 137 void GoToLast() const ; 165 138 bool IsEmpty() const ; 166 139 bool IsEnd() const ; … … 257 230 258 231 /// Count and change present atoms' coordination. 232 void CountAtoms(); 259 233 void CountElements(); 260 234 void CalculateOrbitals(class config &configuration); … … 325 299 bool StoreForcesFile(MoleculeListClass *BondFragments, char *path, int *SortIndex); 326 300 bool CreateMappingLabelsToConfigSequence(int *&SortIndex); 327 bool CreateFatherLookupTable(atom **&LookupTable, int count = 0);328 301 void BreadthFirstSearchAdd(molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem); 329 302 /// -# BOSSANOVA … … 354 327 private: 355 328 int last_atom; //!< number given to last atom 356 mutable internal_iteratorInternalPointer; //!< internal pointer for PointCloud329 mutable atom *InternalPointer; //!< internal pointer for PointCloud 357 330 }; 358 331 -
src/molecule_dynamics.cpp
ra7b761b r8f215d 28 28 gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM); 29 29 gsl_vector *x = gsl_vector_alloc(NDIM); 30 atom * Runner = mol->start; 30 31 atom *Sprinter = NULL; 31 32 Vector trajectory1, trajectory2, normal, TestVector; 32 33 double Norm1, Norm2, tmp, result = 0.; 33 34 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++) 36 38 break; 37 39 // determine normalized trajectories direction vector (n1, n2) … … 40 42 trajectory1.Normalize(); 41 43 Norm1 = trajectory1.Norm(); 42 Sprinter = Params.PermutationMap[ (*iter)->nr]; // find second target point43 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); 44 46 trajectory2.Normalize(); 45 47 Norm2 = trajectory1.Norm(); 46 48 // check whether either is zero() 47 49 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)); 49 51 } else if (Norm1 < MYEPSILON) { 50 52 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); 52 54 trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything 53 55 trajectory1 -= trajectory2; // project the part in norm direction away 54 56 tmp = trajectory1.Norm(); // remaining norm is distance 55 57 } else if (Norm2 < MYEPSILON) { 56 Sprinter = Params.PermutationMap[ (*iter)->nr]; // find second target point58 Sprinter = Params.PermutationMap[Runner->nr]; // find second target point 57 59 trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep); // copy second offset 58 60 trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything … … 64 66 // Log() << Verbose(0) << " and "; 65 67 // 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)); 67 69 // Log() << Verbose(0) << " with distance " << tmp << "." << endl; 68 70 } 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 "; 70 72 // Log() << Verbose(0) << endl; 71 73 // Log() << Verbose(0) << "First Trajectory: "; … … 83 85 gsl_matrix_set(A, 1, i, trajectory2[i]); 84 86 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])); 86 88 } 87 89 // solve the linear system by Householder transformations … … 94 96 trajectory2.Scale(gsl_vector_get(x,1)); 95 97 normal.Scale(gsl_vector_get(x,2)); 96 TestVector = (*iter)->Trajectory.R.at(Params.startstep) + trajectory2 + normal98 TestVector = Runner->Trajectory.R.at(Params.startstep) + trajectory2 + normal 97 99 - (Walker->Trajectory.R.at(Params.startstep) + trajectory1); 98 100 if (TestVector.Norm() < MYEPSILON) { … … 123 125 { 124 126 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)) { 127 131 // 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 "; 129 133 // Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep); 130 134 // Log() << Verbose(0) << ", penalting." << endl; … … 157 161 // go through every atom 158 162 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; 160 166 // first term: distance to target 161 Runner = Params.PermutationMap[ (*iter)->nr]; // find target point162 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))); 163 169 tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem; 164 170 result += Params.PenaltyConstants[0] * tmp; … … 166 172 167 173 // second term: sum of distances to other trajectories 168 result += SumDistanceOfTrajectories( (*iter), this, Params);174 result += SumDistanceOfTrajectories(Walker, this, Params); 169 175 170 176 // third term: penalty for equal targets 171 result += PenalizeEqualTargets( (*iter), this, Params);177 result += PenalizeEqualTargets(Walker, this, Params); 172 178 } 173 179 … … 207 213 void FillDistanceList(molecule *mol, struct EvaluatePotential &Params) 208 214 { 209 for (int i=mol-> getAtomCount(); i--;) {215 for (int i=mol->AtomCount; i--;) { 210 216 Params.DistanceList[i] = new DistanceMap; // is the distance sorted target list per atom 211 217 Params.DistanceList[i]->clear(); 212 218 } 213 219 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) ); 217 228 } 218 229 } … … 226 237 void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params) 227 238 { 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); 234 247 } 235 248 }; … … 272 285 void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params) 273 286 { 274 molecule::const_iterator iter = mol->begin();287 atom *Walker = mol->start; 275 288 DistanceMap::iterator NewBase; 276 289 double Potential = fabs(mol->ConstrainedPotential(Params)); 277 290 278 if (mol->empty()) {279 eLog() << Verbose(1) << "Molecule is empty." << endl;280 return;281 }282 291 while ((Potential) > Params.PenaltyConstants[2]) { 283 PrintPermutationMap(mol-> getAtomCount(), Params);284 iter++;285 if ( iter == mol->end()) // round-robin at the end286 iter = mol->begin();287 if (Params.DoubleList[Params.DistanceIterators[ (*iter)->nr]->second->nr] <= 1) // no need to make those injective that aren't292 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 288 297 continue; 289 298 // 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 <=1299 Potential = TryNextNearestNeighbourForInjectivePermutation(mol, Walker, Potential, Params); 300 } 301 for (int i=mol->AtomCount; i--;) // now each single entry in the DoubleList should be <=1 293 302 if (Params.DoubleList[i] > 1) { 294 303 DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl); … … 329 338 double Potential, OldPotential, OlderPotential; 330 339 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"); 336 345 int round; 337 atom * Sprinter = NULL;346 atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL; 338 347 DistanceMap::iterator Rider, Strider; 339 348 … … 362 371 DoLog(2) && (Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl); 363 372 OlderPotential = OldPotential; 364 molecule::const_iterator iter;365 373 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(); 374 383 break; 375 384 } 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; 377 386 // 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; 382 391 break; 383 392 } 393 Runner = Runner->next; 384 394 } 385 if ( runner != end()) { // we found the other source395 if (Runner != end) { // we found the other source 386 396 // 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++) 389 399 if (Rider->second == Sprinter) 390 400 break; 391 if (Rider != Params.DistanceList[ (*runner)->nr]->end()) { // if we have found one392 //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; 393 403 // exchange both 394 Params.PermutationMap[ (*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // put next farther distance into PermutationMap395 Params.PermutationMap[ (*runner)->nr] = Sprinter; // and hand the old target to its respective owner396 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); 397 407 // calculate the new potential 398 408 //Log() << Verbose(2) << "Checking new potential ..." << endl; … … 400 410 if (Potential > OldPotential) { // we made everything worse! Undo ... 401 411 //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; 403 413 // 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; 405 415 // Undo for Walker 406 Params.DistanceIterators[ (*iter)->nr] = Strider; // take next farther distance target407 //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; 409 419 } else { 410 Params.DistanceIterators[ (*runner)->nr] = Rider; // if successful also move the pointer in the iterator list420 Params.DistanceIterators[Runner->nr] = Rider; // if successful also move the pointer in the iterator list 411 421 DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl); 412 422 OldPotential = Potential; … … 418 428 //Log() << Verbose(0) << endl; 419 429 } 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); 421 431 exit(255); 422 432 } 423 433 } 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! 425 435 } 426 Params.StepList[ (*iter)->nr]++; // take next farther distance target436 Params.StepList[Walker->nr]++; // take next farther distance target 427 437 } 428 } while ( ++iter != end());438 } while (Walker->next != end); 429 439 } while ((OlderPotential - OldPotential) > 1e-3); 430 440 DoLog(1) && (Log() << Verbose(1) << "done." << endl); … … 432 442 433 443 /// free memory and return with evaluated potential 434 for (int i= getAtomCount(); i--;)444 for (int i=AtomCount; i--;) 435 445 Params.DistanceList[i]->clear(); 436 446 Free(&Params.DistanceList); … … 473 483 // Get the Permutation Map by MinimiseConstrainedPotential 474 484 atom **PermutationMap = NULL; 475 atom * Sprinter = NULL;485 atom *Walker = NULL, *Sprinter = NULL; 476 486 if (!MapByIdentity) 477 487 MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem()); 478 488 else { 479 PermutationMap = Malloc<atom *>( getAtomCount(), "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");489 PermutationMap = Malloc<atom *>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap"); 480 490 SetIndexedArrayForEachAtomTo( PermutationMap, &atom::nr ); 481 491 } … … 492 502 mol = World::getInstance().createMolecule(); 493 503 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; 495 507 // add to molecule list 496 Sprinter = mol->AddCopyAtom( (*iter));508 Sprinter = mol->AddCopyAtom(Walker); 497 509 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); 499 511 // add to Trajectories 500 512 //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl; 501 513 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.; 505 517 } 506 518 } … … 510 522 511 523 // 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--; ) 514 526 SortIndex[i] = i; 515 527 status = MoleculePerStep->OutputConfigForListOfFragments(&configuration, SortIndex); … … 555 567 return false; 556 568 } 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); 559 571 performCriticalExit(); 560 572 return false; … … 562 574 // correct Forces 563 575 Velocity.Zero(); 564 for(int i=0;i< getAtomCount();i++)576 for(int i=0;i<AtomCount;i++) 565 577 for(int d=0;d<NDIM;d++) { 566 578 Velocity[d] += Force.Matrix[0][i][d+5]; 567 579 } 568 for(int i=0;i< getAtomCount();i++)580 for(int i=0;i<AtomCount;i++) 569 581 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; 571 583 } 572 584 // solve a constrained potential if we are meant to … … 671 683 delta_alpha = 0.; 672 684 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); 674 686 configuration.alpha += delta_alpha*configuration.Deltat; 675 687 DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl); -
src/molecule_fragmentation.cpp
ra7b761b r8f215d 39 39 int FragmentCount; 40 40 // 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; 43 45 } 44 46 FragmentCount = NoNonHydrogen*(1 << (c*order)); … … 357 359 map<double, pair<int,int> > * ReMapAdaptiveCriteriaListToValue(map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol) 358 360 { 359 atom *Walker = NULL;361 atom *Walker = mol->start; 360 362 map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ; 361 363 DoLog(1) && (Log() << Verbose(1) << "Root candidate list is: " << endl); … … 389 391 bool MarkUpdateCandidates(bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol) 390 392 { 391 atom *Walker = NULL;393 atom *Walker = mol->start; 392 394 int No = -1; 393 395 bool status = false; … … 433 435 bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path) 434 436 { 437 atom *Walker = start; 435 438 bool status = false; 436 439 437 440 // initialize mask list 438 for(int i= getAtomCount();i--;)441 for(int i=AtomCount;i--;) 439 442 AtomMask[i] = false; 440 443 441 444 if (Order < 0) { // adaptive increase of BondOrder per site 442 if (AtomMask[ getAtomCount()] == true) // break after one step445 if (AtomMask[AtomCount] == true) // break after one step 443 446 return false; 444 447 … … 454 457 if (AdaptiveCriteriaList->empty()) { 455 458 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; 457 461 #ifdef ADDHYDROGEN 458 if ( (*iter)->type->Z != 1) // skip hydrogen462 if (Walker->type->Z != 1) // skip hydrogen 459 463 #endif 460 464 { 461 AtomMask[ (*iter)->nr] = true; // include all (non-hydrogen) atoms465 AtomMask[Walker->nr] = true; // include all (non-hydrogen) atoms 462 466 status = true; 463 467 } … … 474 478 Free(&FinalRootCandidates); 475 479 } 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; 477 482 #ifdef ADDHYDROGEN 478 if ( (*iter)->type->Z != 1) // skip hydrogen483 if (Walker->type->Z != 1) // skip hydrogen 479 484 #endif 480 485 { 481 AtomMask[ (*iter)->nr] = true; // include all (non-hydrogen) atoms482 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])) 483 488 status = true; 484 489 } 485 490 } 486 if (( !Order) && (!AtomMask[getAtomCount()])) // single stepping, just check491 if ((Order == 0) && (AtomMask[AtomCount] == false)) // single stepping, just check 487 492 status = true; 488 493 … … 495 500 } 496 501 497 PrintAtomMask(AtomMask, getAtomCount()); // for debugging502 PrintAtomMask(AtomMask, AtomCount); // for debugging 498 503 499 504 return status; … … 511 516 return false; 512 517 } 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--;) 515 520 SortIndex[i] = -1; 516 521 … … 519 524 520 525 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 - failure531 */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 them543 if (count == 0) {544 for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron545 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 fill554 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;574 526 }; 575 527 … … 596 548 MoleculeListClass *BondFragments = NULL; 597 549 int *SortIndex = NULL; 598 int *MinimumRingSize = new int[ getAtomCount()];550 int *MinimumRingSize = new int[AtomCount]; 599 551 int FragmentCounter; 600 552 MoleculeLeafClass *MolecularWalker = NULL; … … 624 576 625 577 // create lookup table for Atom::nr 626 FragmentationToDo = FragmentationToDo && CreateFatherLookupTable( ListOfAtoms, getAtomCount());578 FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(start, end, ListOfAtoms, AtomCount); 627 579 628 580 // === compare it with adjacency file === … … 634 586 635 587 // 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; 638 590 MolecularWalker = Subgraphs; 639 591 FragmentCounter = 0; … … 672 624 // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle 673 625 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; 676 628 FragmentationToDo = false; // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards 677 629 while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) { 678 630 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() 680 632 // ===== 6b. fill RootStack for each subgraph (second adaptivity check) ===== 681 633 Subgraphs->next->FillRootStackForSubgraphs(RootStack, AtomMask, (FragmentCounter = 0)); … … 816 768 bool molecule::ParseOrderAtSiteFromFile(char *path) 817 769 { 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"); 820 772 bool status; 821 773 int AtomNr, value; … … 920 872 atom *OtherFather = NULL; 921 873 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; 928 879 LonelyFlag = true; 929 FatherOfRunner = (*iter)->father; 930 ASSERT(FatherOfRunner,"Atom without father found"); 880 FatherOfRunner = Runner->father; 931 881 if (SonList[FatherOfRunner->nr] != NULL) { // check if this, our father, is present in list 932 882 // create all bonds … … 939 889 // Log() << Verbose(3) << "Adding Bond: "; 940 890 // Log() << Verbose(0) << 941 Leaf->AddBond( (*iter), SonList[OtherFather->nr], (*BondRunner)->BondDegree);891 Leaf->AddBond(Runner, SonList[OtherFather->nr], (*BondRunner)->BondDegree); 942 892 // Log() << Verbose(0) << "." << endl; 943 //NumBonds[ (*iter)->nr]++;893 //NumBonds[Runner->nr]++; 944 894 } else { 945 895 // Log() << Verbose(3) << "Not adding bond, labels in wrong order." << endl; … … 949 899 // Log() << Verbose(0) << ", who has no son in this fragment molecule." << endl; 950 900 #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)) 953 903 exit(1); 954 904 #endif 955 //NumBonds[ (*iter)->nr] += Binder->BondDegree;905 //NumBonds[Runner->nr] += Binder->BondDegree; 956 906 } 957 907 } 958 908 } 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 } 965 914 #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; 969 917 #endif 970 918 } … … 981 929 molecule * molecule::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem) 982 930 { 983 atom **SonList = Calloc<atom*>( getAtomCount(), "molecule::StoreFragmentFromStack: **SonList");931 atom **SonList = Calloc<atom*>(AtomCount, "molecule::StoreFragmentFromStack: **SonList"); 984 932 molecule *Leaf = World::getInstance().createMolecule(); 985 933 … … 1608 1556 FragmentSearch.FragmentSet = new KeySet; 1609 1557 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--;) { 1612 1560 FragmentSearch.ShortestPathList[i] = -1; 1613 1561 } 1614 1562 1615 1563 // Construct the complete KeySet which we need for topmost level only (but for all Roots) 1564 atom *Walker = start; 1616 1565 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); 1619 1569 } 1620 1570 … … 1627 1577 RootKeyNr = RootStack.front(); 1628 1578 RootStack.pop_front(); 1629 atom *Walker = FindAtom(RootKeyNr);1579 Walker = FindAtom(RootKeyNr); 1630 1580 // check cyclic lengths 1631 1581 //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) { … … 1714 1664 Vector Translationvector; 1715 1665 //class StackClass<atom *> *CompStack = NULL; 1716 class StackClass<atom *> *AtomStack = new StackClass<atom *>( getAtomCount());1666 class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount); 1717 1667 bool flag = true; 1718 1668 1719 1669 DoLog(2) && (Log() << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl); 1720 1670 1721 ColorList = Calloc<enum Shading>( getAtomCount(), "molecule::ScanForPeriodicCorrection: *ColorList");1671 ColorList = Calloc<enum Shading>(AtomCount, "molecule::ScanForPeriodicCorrection: *ColorList"); 1722 1672 while (flag) { 1723 1673 // remove bonds that are beyond bonddistance … … 1751 1701 Log() << Verbose(0) << Translationvector << endl; 1752 1702 // apply to all atoms of first component via BFS 1753 for (int i= getAtomCount();i--;)1703 for (int i=AtomCount;i--;) 1754 1704 ColorList[i] = white; 1755 1705 AtomStack->Push(Binder->leftatom); -
src/molecule_geometry.cpp
ra7b761b r8f215d 69 69 70 70 // Log() << Verbose(3) << "Begin of CenterEdge." << endl; 71 molecule::const_iterator iter = begin(); // start at first in list72 if ( iter != end()) {//list not empty?71 atom *ptr = start->next; // start at first in list 72 if (ptr != end) { //list not empty? 73 73 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); 79 80 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); 82 83 } 83 84 } … … 103 104 { 104 105 int Num = 0; 105 molecule::const_iterator iter = begin(); // start at first in list106 atom *ptr = start; // start at first in list 106 107 107 108 Center.Zero(); 108 109 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; 111 113 Num++; 112 Center += (*iter)->x;114 Center += ptr->x; 113 115 } 114 116 Center.Scale(-1./Num); // divide through total number (and sign for direction) … … 123 125 Vector * molecule::DetermineCenterOfAll() const 124 126 { 125 molecule::const_iterator iter = begin(); // start at first in list127 atom *ptr = start->next; // start at first in list 126 128 Vector *a = new Vector(); 127 129 Vector tmp; … … 130 132 a->Zero(); 131 133 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; 134 137 Num += 1.; 135 tmp = (*iter)->x;138 tmp = ptr->x; 136 139 (*a) += tmp; 137 140 } … … 147 150 Vector * molecule::DetermineCenterOfGravity() 148 151 { 149 molecule::const_iterator iter = begin(); // start at first in list152 atom *ptr = start->next; // start at first in list 150 153 Vector *a = new Vector(); 151 154 Vector tmp; … … 154 157 a->Zero(); 155 158 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; 160 164 (*a) += tmp; 161 165 } … … 198 202 void molecule::Scale(const double ** const factor) 199 203 { 200 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 204 atom *ptr = start; 205 206 while (ptr->next != end) { 207 ptr = ptr->next; 201 208 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); 204 211 } 205 212 }; … … 210 217 void molecule::Translate(const Vector *trans) 211 218 { 212 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 219 atom *ptr = start; 220 221 while (ptr->next != end) { 222 ptr = ptr->next; 213 223 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); 216 226 } 217 227 }; … … 249 259 void molecule::DeterminePeriodicCenter(Vector ¢er) 250 260 { 261 atom *Walker = start; 251 262 double * const cell_size = World::getInstance().getDomain(); 252 263 double *matrix = ReturnFullMatrixforSymmetric(cell_size); … … 259 270 Center.Zero(); 260 271 flag = true; 261 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 272 while (Walker->next != end) { 273 Walker = Walker->next; 262 274 #ifdef ADDHYDROGEN 263 if ( (*iter)->type->Z != 1) {275 if (Walker->type->Z != 1) { 264 276 #endif 265 Testvector = (*iter)->x;277 Testvector = Walker->x; 266 278 Testvector.MatrixMultiplication(inversematrix); 267 279 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 nothing280 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 270 282 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]; 272 284 if ((fabs(tmp)) > BondDistance) { 273 285 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); 275 287 if (tmp > 0) 276 288 Translationvector[j] -= 1.; … … 286 298 #ifdef ADDHYDROGEN 287 299 // 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; 291 303 Testvector.MatrixMultiplication(inversematrix); 292 304 Testvector += Translationvector; … … 303 315 Free(&inversematrix); 304 316 305 Center.Scale(1./ static_cast<double>(getAtomCount()));317 Center.Scale(1./(double)AtomCount); 306 318 }; 307 319 … … 313 325 void molecule::PrincipalAxisSystem(bool DoRotate) 314 326 { 327 atom *ptr = start; // start at first in list 315 328 double InertiaTensor[NDIM*NDIM]; 316 329 Vector *CenterOfGravity = DetermineCenterOfGravity(); … … 323 336 324 337 // 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; 327 341 //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]); 337 351 } 338 352 // print InertiaTensor for debugging … … 372 386 373 387 // 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]); 385 402 } 386 403 // print InertiaTensor for debugging … … 406 423 void molecule::Align(Vector *n) 407 424 { 425 atom *ptr = start; 408 426 double alpha, tmp; 409 427 Vector z_axis; … … 416 434 alpha = atan(-n->at(0)/n->at(2)); 417 435 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]; 422 441 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]; 426 445 } 427 446 } … … 433 452 434 453 // rotate on z-y plane 454 ptr = start; 435 455 alpha = atan(-n->at(1)/n->at(2)); 436 456 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]; 441 462 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]; 445 466 } 446 467 } … … 466 487 Vector a,b,c,d; 467 488 struct lsq_params *par = (struct lsq_params *)params; 489 atom *ptr = par->mol->start; 468 490 469 491 // initialize vectors … … 475 497 b[2] = gsl_vector_get(x,5); 476 498 // 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; 480 503 t = c.ScalarProduct(b); // get direction parameter 481 504 d = t*b; // and create vector -
src/molecule_graph.cpp
ra7b761b r8f215d 19 19 #include "World.hpp" 20 20 #include "Helpers/fast_functions.hpp" 21 #include "Helpers/Assert.hpp"22 23 21 24 22 struct BFSAccounting … … 76 74 flip(atom1, atom2); 77 75 Walker = FindAtom(atom1); 78 ASSERT(Walker,"Could not find an atom with the ID given in dbond file");79 76 OtherWalker = FindAtom(atom2); 80 ASSERT(OtherWalker,"Could not find an atom with the ID given in dbond file");81 77 AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices. 82 78 } … … 107 103 atom *Walker = NULL; 108 104 atom *OtherWalker = NULL; 105 atom **AtomMap = NULL; 109 106 int n[NDIM]; 110 107 double MinDistance, MaxDistance; … … 131 128 132 129 // 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.)) { 136 134 DoLog(2) && (Log() << Verbose(2) << "Creating Linked Cell structure ... " << endl); 137 135 LC = new LinkedCell(this, bonddistance); … … 139 137 // create a list to map Tesselpoint::nr to atom * 140 138 DoLog(2) && (Log() << Verbose(2) << "Creating TesselPoint to atom map ... " << endl); 141 142 // set numbers for atoms that can later be used143 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; 146 144 } 147 145 … … 155 153 if (List != NULL) { 156 154 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; 160 157 // 3c. check for possible bond between each atom in this and every one in the 27 cells 161 158 for (n[0] = -1; n[0] <= 1; n[0]++) … … 167 164 for (LinkedCell::LinkedNodes::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) { 168 165 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; 172 170 (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem); 173 const double distance = OtherWalker->x.PeriodicDistanceSquared(Walker->x,cell_size);174 171 const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance); 175 172 // Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl; … … 191 188 } 192 189 } 190 Free(&AtomMap); 193 191 delete (LC); 194 192 DoLog(1) && (Log() << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl); … … 201 199 ActOnAllAtoms( &atom::OutputBondOfAtom ); 202 200 } 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); 204 202 DoLog(0) && (Log() << Verbose(0) << "End of CreateAdjacencyList." << endl); 205 203 if (free_BG) … … 243 241 DoLog(0) && (Log() << Verbose(0) << " done." << endl); 244 242 } 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); 246 244 } 247 245 DoLog(0) && (Log() << Verbose(0) << No << " bonds could not be corrected." << endl); … … 463 461 void DepthFirstSearchAnalysis_Init(struct DFSAccounting &DFS, const molecule * const mol) 464 462 { 465 DFS.AtomStack = new StackClass<atom *> (mol-> getAtomCount());463 DFS.AtomStack = new StackClass<atom *> (mol->AtomCount); 466 464 DFS.CurrentGraphNr = 0; 467 465 DFS.ComponentNumber = 0; … … 504 502 bond *Binder = NULL; 505 503 506 if ( getAtomCount()== 0)504 if (AtomCount == 0) 507 505 return SubGraphs; 508 506 DoLog(0) && (Log() << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl); 509 507 DepthFirstSearchAnalysis_Init(DFS, this); 510 508 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 513 511 // (1) mark all edges unused, empty stack, set atom->GraphNr = -1 for all 514 512 DFS.AtomStack->ClearStack(); … … 550 548 551 549 // 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 on555 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; 556 554 } 557 555 } … … 856 854 if (MinRingSize != -1) { // if rings are present 857 855 // 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 862 861 Walker = Root; 863 862 864 863 //Log() << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl; 865 CyclicStructureAnalysis_BFSToNextCycle(Root, Walker, MinimumRingSize, mol-> getAtomCount());864 CyclicStructureAnalysis_BFSToNextCycle(Root, Walker, MinimumRingSize, mol->AtomCount); 866 865 867 866 } … … 893 892 int MinRingSize = -1; 894 893 895 InitializeBFSAccounting(BFS, getAtomCount());894 InitializeBFSAccounting(BFS, AtomCount); 896 895 897 896 //Log() << Verbose(1) << "Back edge list - "; … … 1135 1134 CurrentBondsOfAtom = -1; // we count one too far due to line end 1136 1135 // parse into structure 1137 if ((AtomNr >= 0) && (AtomNr < getAtomCount())) {1136 if ((AtomNr >= 0) && (AtomNr < AtomCount)) { 1138 1137 Walker = ListOfAtoms[AtomNr]; 1139 1138 while (!line.eof()) … … 1314 1313 AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root); 1315 1314 1316 BreadthFirstSearchAdd_Init(BFS, Root, BondOrder, getAtomCount(), AddedAtomList);1315 BreadthFirstSearchAdd_Init(BFS, Root, BondOrder, AtomCount, AddedAtomList); 1317 1316 1318 1317 // and go on ... Queue always contains all lightgray vertices … … 1370 1369 // fill parent list with sons 1371 1370 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; 1374 1375 // 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 ; 1378 1381 1379 1382 void BuildInducedSubgraph_Finalize(atom **&ParentList) … … 1386 1389 { 1387 1390 bool status = true; 1391 atom *Walker = NULL; 1388 1392 atom *OtherAtom = NULL; 1389 1393 // check each entry of parent list and if ok (one-to-and-onto matching) create bonds 1390 1394 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) { 1394 1400 status = false; 1395 1401 } 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); 1398 1404 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); 1401 1407 } 1402 1408 } … … 1421 1427 bool status = true; 1422 1428 atom **ParentList = NULL; 1429 1423 1430 DoLog(2) && (Log() << Verbose(2) << "Begin of BuildInducedSubgraph." << endl); 1424 BuildInducedSubgraph_Init(ParentList, Father-> getAtomCount());1431 BuildInducedSubgraph_Init(ParentList, Father->AtomCount); 1425 1432 BuildInducedSubgraph_FillParentList(this, Father, ParentList); 1426 1433 status = BuildInducedSubgraph_CreateBondsFromParent(this, Father, ParentList); -
src/molecule_pointcloud.cpp
ra7b761b r8f215d 8 8 #include "atom.hpp" 9 9 #include "config.hpp" 10 #include "info.hpp"11 10 #include "memoryallocator.hpp" 12 11 #include "molecule.hpp" … … 32 31 }; 33 32 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 37 35 */ 38 TesselPoint *molecule::GetPoint() const36 TesselPoint *molecule::GetPoint() const 39 37 { 40 Info FunctionInfo(__func__); 41 return (*InternalPointer); 38 if ((InternalPointer != start) && (InternalPointer != end)) 39 return InternalPointer; 40 else 41 return NULL; 42 42 }; 43 43 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 */ 47 TesselPoint *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 */ 55 int molecule::GetMaxId() const 56 { 57 return last_atom; 58 }; 59 60 /** Go to next atom. 61 * Stops at last one. 46 62 */ 47 63 void molecule::GoToNext() const 48 64 { 49 Info FunctionInfo(__func__); 50 if (InternalPointer != atoms.end()) 51 InternalPointer++; 65 if (InternalPointer != end) 66 InternalPointer = InternalPointer->next; 52 67 }; 53 68 54 /** PointCloud implementation of GoToFirst. 55 * Uses atoms and STL stuff. 69 /** Go to previous atom. 70 * Stops at first one. 71 */ 72 void molecule::GoToPrevious() const 73 { 74 if (InternalPointer->previous != start) 75 InternalPointer = InternalPointer->previous; 76 }; 77 78 /** Goes to first atom. 56 79 */ 57 80 void molecule::GoToFirst() const 58 81 { 59 Info FunctionInfo(__func__); 60 InternalPointer = atoms.begin(); 82 InternalPointer = start->next; 61 83 }; 62 84 63 /** PointCloud implementation of IsEmpty. 64 * Uses atoms and STL stuff. 85 /** Goes to last atom. 86 */ 87 void 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 65 94 */ 66 95 bool molecule::IsEmpty() const 67 96 { 68 Info FunctionInfo(__func__); 69 return (empty()); 97 return (start->next == end); 70 98 }; 71 99 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 74 102 */ 75 103 bool molecule::IsEnd() const 76 104 { 77 Info FunctionInfo(__func__); 78 return (InternalPointer == atoms.end()); 105 return (InternalPointer == end); 79 106 }; 80 81 int molecule::GetMaxId() const {82 return getAtomCount();83 } -
src/molecule_template.hpp
ra7b761b r8f215d 24 24 template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() ) const 25 25 { 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)(); 28 30 } 29 31 }; 30 32 template <typename res> void molecule::ActOnAllVectors( res (Vector::*f)() const ) const 31 33 { 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)(); 34 38 } 35 39 }; … … 37 41 template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T), T t ) const 38 42 { 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); 41 47 } 42 48 }; 43 49 template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T) const, T t ) const 44 50 { 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); 47 55 } 48 56 }; 49 57 template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T&), T &t ) const 50 58 { 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); 53 63 } 54 64 }; 55 65 template <typename res, typename T> void molecule::ActOnAllVectors( res (Vector::*f)(T&) const, T &t ) const 56 66 { 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); 59 71 } 60 72 }; … … 62 74 template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U), T t, U u ) const 63 75 { 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); 66 80 } 67 81 }; 68 82 template <typename res, typename T, typename U> void molecule::ActOnAllVectors( res (Vector::*f)(T, U) const, T t, U u ) const 69 83 { 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); 72 88 } 73 89 }; … … 75 91 template <typename res, typename T, typename U, typename V> void molecule::ActOnAllVectors( res (Vector::*f)(T, U, V), T t, U u, V v) const 76 92 { 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); 79 97 } 80 98 }; 81 99 template <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 82 100 { 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); 85 105 } 86 106 }; … … 92 112 { 93 113 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)(); 96 118 } 97 119 return result; … … 100 122 { 101 123 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)(); 104 128 } 105 129 return result; … … 109 133 { 110 134 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); 113 139 } 114 140 return result; … … 117 143 { 118 144 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); 121 149 } 122 150 return result; … … 129 157 template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *)) const 130 158 { 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); 133 163 } 134 164 }; 135 165 template <typename res> void molecule::ActWithEachAtom( res (molecule::*f)(atom *) const) const 136 166 { 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); 139 171 } 140 172 }; … … 145 177 template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) , molecule *copy) const 146 178 { 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); 149 183 } 150 184 }; 151 185 template <typename res> void molecule::ActOnCopyWithEachAtom( res (molecule::*f)(atom *) const, molecule *copy) const 152 186 { 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); 155 191 } 156 192 }; … … 161 197 template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () ) const 162 198 { 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); 166 204 } 167 205 }; 168 206 template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) () const ) const 169 207 { 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); 173 213 } 174 214 }; 175 215 template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const , molecule *copy, bool (atom::*condition) () ) const 176 216 { 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); 180 222 } 181 223 }; 182 224 template <typename res> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) () const ) const 183 225 { 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); 187 231 } 188 232 }; … … 190 234 template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T), T t ) const 191 235 { 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); 195 241 } 196 242 }; 197 243 template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) , molecule *copy, bool (atom::*condition) (T) const, T t ) const 198 244 { 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); 202 250 } 203 251 }; 204 252 template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T), T t ) const 205 253 { 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); 209 259 } 210 260 }; 211 261 template <typename res, typename T> void molecule::ActOnCopyWithEachAtomIfTrue( res (molecule::*f)(atom *) const, molecule *copy, bool (atom::*condition) (T) const, T t ) const 212 262 { 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); 216 268 } 217 269 }; … … 219 271 template <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 220 272 { 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); 224 278 } 225 279 }; 226 280 template <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 227 281 { 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); 231 287 } 232 288 }; 233 289 template <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 234 290 { 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); 238 296 } 239 297 }; 240 298 template <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 241 299 { 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); 245 305 } 246 306 }; … … 248 308 template <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 249 309 { 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); 253 315 } 254 316 }; 255 317 template <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 256 318 { 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); 260 324 } 261 325 }; 262 326 template <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 263 327 { 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); 267 333 } 268 334 }; 269 335 template <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 270 336 { 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); 274 342 } 275 343 }; … … 280 348 template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)()) const 281 349 { 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)(); 284 354 } 285 355 }; 286 356 template <typename res, typename typ> void molecule::ActOnAllAtoms( res (typ::*f)() const) const 287 357 { 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)(); 290 362 } 291 363 }; … … 293 365 template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T), T t ) const 294 366 { 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); 297 371 } 298 372 }; 299 373 template <typename res, typename typ, typename T> void molecule::ActOnAllAtoms( res (typ::*f)(T) const, T t ) const 300 374 { 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); 303 379 } 304 380 }; … … 306 382 template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U), T t, U u ) const 307 383 { 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); 310 388 } 311 389 }; 312 390 template <typename res, typename typ, typename T, typename U> void molecule::ActOnAllAtoms( res (typ::*f)(T, U) const, T t, U u ) const 313 391 { 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); 316 396 } 317 397 }; … … 319 399 template <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 320 400 { 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); 323 405 } 324 406 }; 325 407 template <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 326 408 { 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); 329 413 } 330 414 }; … … 332 416 template <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 333 417 { 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); 336 422 } 337 423 }; 338 424 template <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 339 425 { 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); 342 430 } 343 431 }; … … 348 436 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *) ) const 349 437 { 438 atom *Walker = start; 350 439 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); 353 443 } 354 444 }; 355 445 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T value ) const 356 446 { 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); 359 451 } 360 452 }; 361 453 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, void (*Setor)(T *, T *), T *value ) const 362 454 { 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); 365 459 } 366 460 }; … … 368 462 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *) ) const 369 463 { 464 atom *Walker = start; 370 465 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); 373 469 } 374 470 }; 375 471 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T value ) const 376 472 { 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); 379 477 } 380 478 }; 381 479 template <typename T> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int element::*index, void (*Setor)(T *, T *), T *value ) const 382 480 { 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); 385 485 } 386 486 }; … … 388 488 template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ atom::*value ) const 389 489 { 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); 392 494 } 393 495 }; 394 496 template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ atom::*value ) const 395 497 { 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); 398 502 } 399 503 }; 400 504 template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &), typ &vect ) const 401 505 { 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); 404 510 } 405 511 }; 406 512 template <typename T, typename typ> void molecule::SetIndexedArrayForEachAtomTo ( T *array, int ParticleInfo::*index, T (atom::*Setor)(typ &) const, typ &vect ) const 407 513 { 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); 410 518 } 411 519 }; 412 520 template <typename T, typename typ, typename typ2> void molecule::SetAtomValueToIndexedArray ( T *array, int typ::*index, T typ2::*value ) const 413 521 { 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; 417 527 } 418 528 }; … … 420 530 template <typename T, typename typ> void molecule::SetAtomValueToValue ( T value, T typ::*ptr ) const 421 531 { 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; 425 537 } 426 538 }; -
src/moleculelist.cpp
ra7b761b r8f215d 20 20 #include "memoryallocator.hpp" 21 21 #include "periodentafel.hpp" 22 #include " Helpers/Assert.hpp"22 #include "World.hpp" 23 23 24 24 /*********************************** Functions for class MoleculeListClass *************************/ … … 69 69 int Count, Counter, aCounter, bCounter; 70 70 int flag; 71 atom *aWalker = NULL; 72 atom *bWalker = NULL; 71 73 72 74 // sort each atom list and put the numbers into a list, then go through 73 75 //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) { 80 77 return -1; 81 78 } else { 82 if ( mol1->getAtomCount() > mol2->getAtomCount())79 if ((**(molecule **) a).AtomCount > (**(molecule **) b).AtomCount) 83 80 return +1; 84 81 else { 85 Count = mol1->getAtomCount();82 Count = (**(molecule **) a).AtomCount; 86 83 aList = new int[Count]; 87 84 bList = new int[Count]; 88 85 89 86 // fill the lists 87 aWalker = (**(molecule **) a).start; 88 bWalker = (**(molecule **) b).start; 90 89 Counter = 0; 91 90 aCounter = 0; 92 91 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) 98 96 aList[Counter] = Count + (aCounter++); 99 97 else 100 aList[Counter] = (*aiter)->GetTrueFather()->nr;101 if ( (*biter)->GetTrueFather() == NULL)98 aList[Counter] = aWalker->GetTrueFather()->nr; 99 if (bWalker->GetTrueFather() == NULL) 102 100 bList[Counter] = Count + (bCounter++); 103 101 else 104 bList[Counter] = (*biter)->GetTrueFather()->nr;102 bList[Counter] = bWalker->GetTrueFather()->nr; 105 103 Counter++; 106 104 } 107 105 // check if AtomCount was for real 108 106 flag = 0; 109 if ((a iter == mol1->end()) && (biter != mol2->end())) {107 if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) { 110 108 flag = -1; 111 109 } else { 112 if ((a iter != mol1->end()) && (biter == mol2->end()))110 if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **) b).end)) 113 111 flag = 1; 114 112 } … … 144 142 void MoleculeListClass::Enumerate(ostream *out) 145 143 { 144 atom *Walker = NULL; 146 145 periodentafel *periode = World::getInstance().getPeriode(); 147 146 std::map<atomicNumber_t,unsigned int> counts; … … 159 158 // count atoms per element and determine size of bounding sphere 160 159 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); 165 166 } 166 167 // 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"; 168 169 169 170 std::map<atomicNumber_t,unsigned int>::reverse_iterator iter; … … 201 202 202 203 // 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); 208 211 } 209 212 … … 225 228 226 229 // 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); 230 236 Walker->father = Walker; 231 237 } … … 324 330 325 331 // 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++) 328 335 CopyAtoms[i] = NULL; 329 336 330 337 // for each of the source atoms check whether we are in- or outside and add copy atom 338 atom *Walker = srcmol->start; 331 339 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]); 337 346 nr++; 338 347 } else { … … 340 349 } 341 350 } 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."); 343 352 344 353 // go through all bonds and add as well … … 373 382 bool MoleculeListClass::AddHydrogenCorrection(char *path) 374 383 { 384 atom *Walker = NULL; 385 atom *Runner = NULL; 375 386 bond *Binder = NULL; 376 387 double ***FitConstant = NULL, **correction = NULL; … … 477 488 correction[k][j] = 0.; 478 489 // 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; 485 500 // 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!) 488 503 // 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; 491 506 for (int k = 0; k < a; k++) { 492 507 for (int j = 0; j < b; j++) { … … 562 577 ofstream ForcesFile; 563 578 stringstream line; 579 atom *Walker = NULL; 564 580 periodentafel *periode=World::getInstance().getPeriode(); 565 581 … … 574 590 periodentafel::const_iterator elemIter; 575 591 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 580 598 //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"; 582 600 } else 583 601 // otherwise a -1 to indicate an added saturation hydrogen … … 615 633 bool result = true; 616 634 bool intermediateResult = true; 635 atom *Walker = NULL; 617 636 Vector BoxDimension; 618 637 char *FragmentNumber = NULL; … … 655 674 // list atoms in fragment for debugging 656 675 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() << " "); 659 680 } 660 681 DoLog(0) && (Log() << Verbose(0) << endl); … … 736 757 { 737 758 molecule *mol = World::getInstance().createMolecule(); 759 atom *Walker = NULL; 760 atom *Advancer = NULL; 738 761 bond *Binder = NULL; 739 762 bond *Stepper = NULL; … … 741 764 for (MoleculeList::iterator MolRunner = ListOfMolecules.begin(); !ListOfMolecules.empty(); MolRunner = ListOfMolecules.begin()) { 742 765 // 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 748 774 } 749 775 // remove all bonds … … 799 825 // 4b. create and fill map of which atom is associated to which connected molecule (note, counting starts at 1) 800 826 int FragmentCounter = 0; 801 int *MolMap = Calloc<int>(mol-> getAtomCount(), "config::Load() - *MolMap");827 int *MolMap = Calloc<int>(mol->AtomCount, "config::Load() - *MolMap"); 802 828 MoleculeLeafClass *MolecularWalker = Subgraphs; 829 Walker = NULL; 803 830 while (MolecularWalker->next != NULL) { 804 831 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; 807 836 } 808 837 FragmentCounter++; … … 810 839 811 840 // 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); 815 846 performCriticalExit(); 816 847 } 817 FragmentCounter = MolMap[ (*iter)->nr];848 FragmentCounter = MolMap[Walker->nr]; 818 849 if (FragmentCounter != 0) { 819 DoLog(3) && (Log() << Verbose(3) << "Re-linking " << * *iter << "..." << endl);820 molecules[FragmentCounter-1]->AddAtom((*iter)); // counting starts at 1821 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 822 853 } 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); 824 855 performCriticalExit(); 825 856 } … … 829 860 while (mol->first->next != mol->last) { 830 861 Binder = mol->first->next; 831 const atom * constWalker = Binder->leftatom;862 Walker = Binder->leftatom; 832 863 unlink(Binder); 833 864 link(Binder,molecules[MolMap[Walker->nr]-1]->last); // counting starts at 1 … … 851 882 int MoleculeListClass::CountAllAtoms() const 852 883 { 884 atom *Walker = NULL; 853 885 int AtomNo = 0; 854 886 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 } 856 892 } 857 893 return AtomNo; … … 1028 1064 bool MoleculeLeafClass::FillBondStructureFromReference(const molecule * const reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList) 1029 1065 { 1066 atom *Walker = NULL; 1030 1067 atom *OtherWalker = NULL; 1031 1068 atom *Father = NULL; … … 1035 1072 DoLog(1) && (Log() << Verbose(1) << "Begin of FillBondStructureFromReference." << endl); 1036 1073 // fill ListOfLocalAtoms if NULL was given 1037 if (!FillListOfLocalAtoms(ListOfLocalAtoms, FragmentCounter, reference-> getAtomCount(), FreeList)) {1074 if (!FillListOfLocalAtoms(ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) { 1038 1075 DoLog(1) && (Log() << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl); 1039 1076 return false; … … 1051 1088 } 1052 1089 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(); 1055 1094 AtomNo = Father->nr; // global id of the current walker 1056 1095 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 walker1096 OtherWalker = ListOfLocalAtoms[FragmentCounter][(*Runner)->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker 1058 1097 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); 1061 1100 } 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); 1063 1102 status = false; 1064 1103 } … … 1087 1126 bool MoleculeLeafClass::FillRootStackForSubgraphs(KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter) 1088 1127 { 1089 atom * Father = NULL;1128 atom *Walker = NULL, *Father = NULL; 1090 1129 1091 1130 if (RootStack != NULL) { … … 1093 1132 if (&(RootStack[FragmentCounter]) != NULL) { 1094 1133 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(); 1097 1138 if (AtomMask[Father->nr]) // apply mask 1098 1139 #ifdef ADDHYDROGEN 1099 if ( (*iter)->type->Z != 1) // skip hydrogen1140 if (Walker->type->Z != 1) // skip hydrogen 1100 1141 #endif 1101 RootStack[FragmentCounter].push_front( (*iter)->nr);1142 RootStack[FragmentCounter].push_front(Walker->nr); 1102 1143 } 1103 1144 if (next != NULL) … … 1138 1179 1139 1180 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); 1141 1182 FreeList = FreeList && true; 1142 1183 } … … 1162 1203 DoLog(1) && (Log() << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl); 1163 1204 // fill ListOfLocalAtoms if NULL was given 1164 if (!FillListOfLocalAtoms(ListOfLocalAtoms, FragmentCounter, reference-> getAtomCount(), FreeList)) {1205 if (!FillListOfLocalAtoms(ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) { 1165 1206 DoLog(1) && (Log() << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl); 1166 1207 return false; -
src/tesselation.cpp
ra7b761b r8f215d 21 21 #include "Plane.hpp" 22 22 #include "Exceptions/LinearDependenceException.hpp" 23 #include "Helpers/Assert.hpp"24 25 23 #include "Helpers/Assert.hpp" 26 24 … … 360 358 // get ascending order of endpoints 361 359 PointMap OrderMap; 362 for (int i = 0; i < 3; i++) {360 for (int i = 0; i < 3; i++) 363 361 // for all three lines 364 362 for (int j = 0; j < 2; j++) { // for both endpoints … … 366 364 // and we don't care whether insertion fails 367 365 } 368 }369 366 // set endpoints 370 367 int Counter = 0; … … 375 372 Counter++; 376 373 } 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 ; 380 380 381 381 /** Destructor of BoundaryTriangleSet. … … 1232 1232 ; 1233 1233 1234 /** PointCloud implementation of GetTerminalPoint. 1235 * Uses PointsOnBoundary and STL stuff. 1236 */ 1237 TesselPoint * 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 1234 1246 /** PointCloud implementation of GoToNext. 1235 1247 * Uses PointsOnBoundary and STL stuff. … … 1243 1255 ; 1244 1256 1257 /** PointCloud implementation of GoToPrevious. 1258 * Uses PointsOnBoundary and STL stuff. 1259 */ 1260 void Tesselation::GoToPrevious() const 1261 { 1262 Info FunctionInfo(__func__); 1263 if (InternalPointer != PointsOnBoundary.begin()) 1264 InternalPointer--; 1265 } 1266 ; 1267 1245 1268 /** PointCloud implementation of GoToFirst. 1246 1269 * Uses PointsOnBoundary and STL stuff. … … 1250 1273 Info FunctionInfo(__func__); 1251 1274 InternalPointer = PointsOnBoundary.begin(); 1275 } 1276 ; 1277 1278 /** PointCloud implementation of GoToLast. 1279 * Uses PointsOnBoundary and STL stuff. 1280 */ 1281 void Tesselation::GoToLast() const 1282 { 1283 Info FunctionInfo(__func__); 1284 InternalPointer = PointsOnBoundary.end(); 1285 InternalPointer--; 1252 1286 } 1253 1287 ; … … 2559 2593 // fill the set of neighbours 2560 2594 TesselPointSet SetOfNeighbours; 2561 2562 2595 SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node); 2563 2596 for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++) -
src/tesselation.hpp
ra7b761b r8f215d 244 244 virtual Vector *GetCenter() const { return NULL; }; 245 245 virtual TesselPoint *GetPoint() const { return NULL; }; 246 virtual TesselPoint *GetTerminalPoint() const { return NULL; }; 246 247 virtual int GetMaxId() const { return 0; }; 247 248 virtual void GoToNext() const {}; 249 virtual void GoToPrevious() const {}; 248 250 virtual void GoToFirst() const {}; 251 virtual void GoToLast() const {}; 249 252 virtual bool IsEmpty() const { return true; }; 250 253 virtual bool IsEnd() const { return true; }; … … 360 363 virtual Vector *GetCenter(ofstream *out) const; 361 364 virtual TesselPoint *GetPoint() const; 365 virtual TesselPoint *GetTerminalPoint() const; 362 366 virtual void GoToNext() const; 367 virtual void GoToPrevious() const; 363 368 virtual void GoToFirst() const; 369 virtual void GoToLast() const; 364 370 virtual bool IsEmpty() const; 365 371 virtual bool IsEnd() const; -
src/tesselationhelpers.cpp
ra7b761b r8f215d 838 838 int i=cloud->GetMaxId(); 839 839 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++) 841 841 LookupList[i] = -1; 842 }843 842 844 843 // print atom coordinates 845 844 int Counter = 1; 846 845 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++) { 848 847 Walker = target->second->node; 849 848 LookupList[Walker->nr] = Counter++; -
src/unittests/AnalysisCorrelationToPointUnitTest.cpp
ra7b761b r8f215d 25 25 #include "periodentafel.hpp" 26 26 #include "tesselation.hpp" 27 #include "World.hpp" 27 28 28 29 #ifdef HAVE_TESTRUNNER … … 79 80 80 81 // check that TestMolecule was correctly constructed 81 CPPUNIT_ASSERT_EQUAL( TestMolecule-> getAtomCount(), 4 );82 CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 ); 82 83 83 84 TestList = World::getInstance().getMolecules(); -
src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
ra7b761b r8f215d 26 26 #include "tesselation.hpp" 27 27 #include "World.hpp" 28 #include "Helpers/Assert.hpp"29 28 30 29 #include "Helpers/Assert.hpp" … … 41 40 void AnalysisCorrelationToSurfaceUnitTest::setUp() 42 41 { 43 ASSERT_DO(Assert::Throw);42 //ASSERT_DO(Assert::Throw); 44 43 45 44 atom *Walker = NULL; … … 72 71 // construct molecule (tetraeder of hydrogens) base 73 72 TestSurfaceMolecule = World::getInstance().createMolecule(); 74 75 73 Walker = World::getInstance().createAtom(); 76 74 Walker->type = hydrogen; 77 75 *Walker->node = Vector(1., 0., 1. ); 78 TestSurfaceMolecule->AddAtom(Walker); 79 76 77 TestSurfaceMolecule->AddAtom(Walker); 80 78 Walker = World::getInstance().createAtom(); 81 79 Walker->type = hydrogen; 82 80 *Walker->node = Vector(0., 1., 1. ); 83 81 TestSurfaceMolecule->AddAtom(Walker); 84 85 82 Walker = World::getInstance().createAtom(); 86 83 Walker->type = hydrogen; 87 84 *Walker->node = Vector(1., 1., 0. ); 88 85 TestSurfaceMolecule->AddAtom(Walker); 89 90 86 Walker = World::getInstance().createAtom(); 91 87 Walker->type = hydrogen; … … 94 90 95 91 // check that TestMolecule was correctly constructed 96 CPPUNIT_ASSERT_EQUAL( TestSurfaceMolecule-> getAtomCount(), 4 );92 CPPUNIT_ASSERT_EQUAL( TestSurfaceMolecule->AtomCount, 4 ); 97 93 98 94 TestList = World::getInstance().getMolecules(); … … 111 107 *Walker->node = Vector(4., 0., 4. ); 112 108 TestSurfaceMolecule->AddAtom(Walker); 113 114 109 Walker = World::getInstance().createAtom(); 115 110 Walker->type = carbon; 116 111 *Walker->node = Vector(0., 4., 4. ); 117 112 TestSurfaceMolecule->AddAtom(Walker); 118 119 113 Walker = World::getInstance().createAtom(); 120 114 Walker->type = carbon; 121 115 *Walker->node = Vector(4., 4., 0. ); 122 116 TestSurfaceMolecule->AddAtom(Walker); 123 124 117 // add inner atoms 125 118 Walker = World::getInstance().createAtom(); … … 127 120 *Walker->node = Vector(0.5, 0.5, 0.5 ); 128 121 TestSurfaceMolecule->AddAtom(Walker); 129 130 122 TestSurfaceMolecule->ActiveFlag = true; 131 123 TestList->insert(TestSurfaceMolecule); … … 158 150 void AnalysisCorrelationToSurfaceUnitTest::SurfaceTest() 159 151 { 160 CPPUNIT_ASSERT_EQUAL( 4, TestSurfaceMolecule-> getAtomCount());152 CPPUNIT_ASSERT_EQUAL( 4, TestSurfaceMolecule->AtomCount ); 161 153 CPPUNIT_ASSERT_EQUAL( (size_t)2, TestList->ListOfMolecules.size() ); 162 154 CPPUNIT_ASSERT_EQUAL( (size_t)4, Surface->PointsOnBoundary.size() ); -
src/unittests/AnalysisPairCorrelationUnitTest.cpp
ra7b761b r8f215d 79 79 80 80 // check that TestMolecule was correctly constructed 81 CPPUNIT_ASSERT_EQUAL( TestMolecule-> getAtomCount(), 4 );81 CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 ); 82 82 83 83 TestList = World::getInstance().getMolecules(); -
src/unittests/CountBondsUnitTest.cpp
ra7b761b r8f215d 102 102 103 103 // 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 ); 106 110 107 111 // create a small file with table -
src/unittests/LinkedCellUnitTest.cpp
ra7b761b r8f215d 69 69 70 70 // 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 ); 72 74 }; 73 75 … … 195 197 { 196 198 // 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) ); 199 203 } 200 204 201 205 // 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); 207 211 208 212 // 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); 214 218 }; 215 219 … … 283 287 size = ListOfPoints->size(); 284 288 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); 288 294 size--; 289 295 CPPUNIT_ASSERT_EQUAL( size, ListOfPoints->size() ); … … 300 306 size=ListOfPoints->size(); 301 307 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); 305 313 size--; 306 314 CPPUNIT_ASSERT_EQUAL( size, ListOfPoints->size() ); … … 318 326 size=ListOfPoints->size(); 319 327 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); 322 332 size--; 323 333 CPPUNIT_ASSERT_EQUAL( size, ListOfPoints->size() ); … … 345 355 size = ListOfPoints->size(); 346 356 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); 350 362 size--; 351 363 CPPUNIT_ASSERT_EQUAL( size, ListOfPoints->size() ); -
src/unittests/ObserverTest.cpp
ra7b761b r8f215d 11 11 #include <cppunit/extensions/TestFactoryRegistry.h> 12 12 #include <cppunit/ui/text/TestRunner.h> 13 #include <set>14 13 15 14 #include "Patterns/Observer.hpp" 16 #include "Patterns/ObservedIterator.hpp"17 15 #include "Helpers/Assert.hpp" 18 16 … … 164 162 bool wasNotified; 165 163 }; 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 204 164 205 165 /******************* actuall tests ***************/ … … 213 173 blockObservable = new BlockObservable(); 214 174 notificationObservable = new NotificationObservable(); 215 collection = new ObservableCollection(5);216 175 217 176 observer1 = new UpdateCountObserver(); … … 222 181 notificationObserver1 = new NotificationObserver(notificationObservable->notification1); 223 182 notificationObserver2 = new NotificationObserver(notificationObservable->notification2); 183 224 184 } 225 185 … … 231 191 delete blockObservable; 232 192 delete notificationObservable; 233 delete collection;234 193 235 194 delete observer1; … … 318 277 blockObservable->changeMethod2(); 319 278 blockObservable->noChangeMethod(); 320 }321 322 void ObserverTest::iteratorTest(){323 int i = 0;324 // test the general iterator methods325 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 of339 // 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 propagated345 CPPUNIT_ASSERT_EQUAL( 0, observer1->updates);346 }347 // After the Iterator has died the propagation should take place348 CPPUNIT_ASSERT_EQUAL( 1, observer1->updates);349 350 // when using a const_iterator no changes should be propagated351 for(ObservableCollection::const_iterator iter = collection->begin(); iter!=collection->end();++iter);352 CPPUNIT_ASSERT_EQUAL( 1, observer1->updates);353 collection->signOff(observer1);354 279 } 355 280 -
src/unittests/ObserverTest.hpp
ra7b761b r8f215d 17 17 class CallObservable; 18 18 class SuperObservable; 19 class ObservableCollection;20 19 class BlockObservable; 21 20 class NotificationObservable; 21 22 22 23 23 class ObserverTest : public CppUnit::TestFixture … … 29 29 CPPUNIT_TEST ( doesNotifyTest ); 30 30 CPPUNIT_TEST ( doesReportTest ); 31 CPPUNIT_TEST ( iteratorTest );32 31 CPPUNIT_TEST ( CircleDetectionTest ); 33 32 CPPUNIT_TEST_SUITE_END(); … … 42 41 void doesNotifyTest(); 43 42 void doesReportTest(); 44 void iteratorTest();45 43 void CircleDetectionTest(); 46 44 … … 60 58 SuperObservable *superObservable; 61 59 NotificationObservable *notificationObservable; 62 ObservableCollection *collection;63 64 60 }; 65 61 -
src/unittests/analysisbondsunittest.cpp
ra7b761b r8f215d 25 25 #include "molecule.hpp" 26 26 #include "periodentafel.hpp" 27 #include "World.hpp"28 27 29 28 #ifdef HAVE_TESTRUNNER … … 90 89 91 90 // check that TestMolecule was correctly constructed 92 CPPUNIT_ASSERT_EQUAL( TestMolecule-> getAtomCount(), 5 );91 CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 5 ); 93 92 94 93 // create a small file with table -
src/unittests/bondgraphunittest.cpp
ra7b761b r8f215d 89 89 90 90 // check that TestMolecule was correctly constructed 91 CPPUNIT_ASSERT_EQUAL( TestMolecule-> getAtomCount(), 4 );91 CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 ); 92 92 93 93 // create a small file with table … … 120 120 }; 121 121 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 130 122 /** UnitTest for BondGraphTest::LoadBondLengthTable(). 131 123 */ … … 142 134 void BondGraphTest::ConstructGraphFromTableTest() 143 135 { 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 ); 147 139 CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(*filename) ); 148 140 CPPUNIT_ASSERT_EQUAL( true , BG->ConstructBondGraph(TestMolecule) ); 149 CPPUNIT_ASSERT_EQUAL( true , (*Walker)->IsBondedTo((*Runner)) );141 CPPUNIT_ASSERT_EQUAL( true , Walker->IsBondedTo(Runner) ); 150 142 }; 151 143 … … 154 146 void BondGraphTest::ConstructGraphFromCovalentRadiiTest() 155 147 { 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 ); 160 151 CPPUNIT_ASSERT_EQUAL( false , BG->LoadBondLengthTable(*dummyname) ); 161 152 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) ); 165 154 }; 166 155 -
src/unittests/bondgraphunittest.hpp
ra7b761b r8f215d 22 22 { 23 23 CPPUNIT_TEST_SUITE( BondGraphTest) ; 24 CPPUNIT_TEST ( SetupTest );25 24 CPPUNIT_TEST ( LoadTableTest ); 26 25 CPPUNIT_TEST ( ConstructGraphFromTableTest ); … … 31 30 void setUp(); 32 31 void tearDown(); 33 void SetupTest();34 32 void LoadTableTest(); 35 33 void ConstructGraphFromTableTest(); -
src/unittests/listofbondsunittest.cpp
ra7b761b r8f215d 74 74 75 75 // check that TestMolecule was correctly constructed 76 CPPUNIT_ASSERT_EQUAL( TestMolecule-> getAtomCount(), 4 );76 CPPUNIT_ASSERT_EQUAL( TestMolecule->AtomCount, 4 ); 77 77 78 78 }; … … 90 90 }; 91 91 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 101 92 /** Unit Test of molecule::AddBond() 102 93 * … … 105 96 { 106 97 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; 111 100 CPPUNIT_ASSERT( atom1 != NULL ); 112 101 CPPUNIT_ASSERT( atom2 != NULL ); … … 135 124 { 136 125 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; 141 128 CPPUNIT_ASSERT( atom1 != NULL ); 142 129 CPPUNIT_ASSERT( atom2 != NULL ); … … 163 150 { 164 151 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; 171 155 CPPUNIT_ASSERT( atom1 != NULL ); 172 156 CPPUNIT_ASSERT( atom2 != NULL ); … … 205 189 { 206 190 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; 211 193 CPPUNIT_ASSERT( atom1 != NULL ); 212 194 CPPUNIT_ASSERT( atom2 != NULL ); … … 233 215 { 234 216 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; 239 219 CPPUNIT_ASSERT( atom1 != NULL ); 240 220 CPPUNIT_ASSERT( atom2 != NULL ); … … 260 240 { 261 241 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; 266 244 CPPUNIT_ASSERT( atom1 != NULL ); 267 245 CPPUNIT_ASSERT( atom2 != NULL ); -
src/unittests/listofbondsunittest.hpp
ra7b761b r8f215d 20 20 { 21 21 CPPUNIT_TEST_SUITE( ListOfBondsTest) ; 22 CPPUNIT_TEST ( SetupTest );23 22 CPPUNIT_TEST ( AddingBondTest ); 24 23 CPPUNIT_TEST ( RemovingBondTest ); … … 32 31 void setUp(); 33 32 void tearDown(); 34 void SetupTest();35 33 void AddingBondTest(); 36 34 void RemovingBondTest();
Note:
See TracChangeset
for help on using the changeset viewer.