Changes in src/molecule.cpp [a67d19:49e1ae]
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/molecule.cpp
ra67d19 r49e1ae 23 23 #include "tesselation.hpp" 24 24 #include "vector.hpp" 25 #include "World.hpp"26 25 27 26 /************************************* Functions for class molecule *********************************/ … … 46 45 for(int i=MAX_ELEMENTS;i--;) 47 46 ElementsInMolecule[i] = 0; 48 strcpy(name,World::get()->DefaultName); 47 cell_size[0] = cell_size[2] = cell_size[5]= 20.; 48 cell_size[1] = cell_size[3] = cell_size[4]= 0.; 49 strcpy(name,"none"); 49 50 }; 50 51 … … 158 159 double *matrix = NULL; 159 160 bond *Binder = NULL; 160 double * const cell_size = World::get()->cell_size;161 161 162 162 // Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl; … … 194 194 BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1]; 195 195 if (BondRescale == -1) { 196 DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);196 eLog() << Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl; 197 197 return false; 198 198 BondRescale = bondlength; … … 237 237 SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin); 238 238 } else { 239 DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name);239 eLog() << Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name; 240 240 } 241 241 } … … 274 274 bondangle = TopOrigin->type->HBondAngle[1]; 275 275 if (bondangle == -1) { 276 DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl);276 eLog() << Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl; 277 277 return false; 278 278 bondangle = 0; … … 396 396 break; 397 397 default: 398 DoeLog(1) && (eLog()<< Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl);398 eLog() << Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl; 399 399 AllWentWell = false; 400 400 break; … … 429 429 input = new istringstream(line); 430 430 *input >> NumberOfAtoms; 431 DoLog(0) && (Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl);431 Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl; 432 432 getline(xyzfile,line,'\n'); // Read comment 433 DoLog(1) && (Log() << Verbose(1) << "Comment: " << line << endl);433 Log() << Verbose(1) << "Comment: " << line << endl; 434 434 435 435 if (MDSteps == 0) // no atoms yet present … … 447 447 Walker->type = elemente->FindElement(shorthand); 448 448 if (Walker->type == NULL) { 449 DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");449 eLog() << Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H."; 450 450 Walker->type = elemente->FindElement(1); 451 451 } … … 543 543 add(Binder, last); 544 544 } else { 545 DoeLog(1) && (eLog()<< Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl);545 eLog() << Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl; 546 546 } 547 547 return Binder; … … 555 555 bool molecule::RemoveBond(bond *pointer) 556 556 { 557 // DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);557 //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl; 558 558 pointer->leftatom->RegisterBond(pointer); 559 559 pointer->rightatom->RegisterBond(pointer); … … 569 569 bool molecule::RemoveBonds(atom *BondPartner) 570 570 { 571 // DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl);571 //eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl; 572 572 BondList::const_iterator ForeRunner; 573 573 while (!BondPartner->ListOfBonds.empty()) { … … 603 603 void molecule::SetBoxDimension(Vector *dim) 604 604 { 605 double * const cell_size = World::get()->cell_size;606 605 cell_size[0] = dim->x[0]; 607 606 cell_size[1] = 0.; … … 622 621 AtomCount--; 623 622 } else 624 DoeLog(1) && (eLog()<< Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl);623 eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl; 625 624 if (ElementsInMolecule[pointer->type->Z] == 0) // was last atom of this element? 626 625 ElementCount--; … … 640 639 ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element 641 640 else 642 DoeLog(1) && (eLog()<< Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl);641 eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl; 643 642 if (ElementsInMolecule[pointer->type->Z] == 0) // was last atom of this element? 644 643 ElementCount--; … … 665 664 return walker; 666 665 } else { 667 DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl);666 Log() << Verbose(0) << "Atom not found in list." << endl; 668 667 return NULL; 669 668 } … … 681 680 //mol->Output((ofstream *)&cout); 682 681 //Log() << Verbose(0) << "===============================================" << endl; 683 DoLog(0) && (Log() << Verbose(0) << text);682 Log() << Verbose(0) << text; 684 683 cin >> No; 685 684 ion = this->FindAtom(No); … … 694 693 bool molecule::CheckBounds(const Vector *x) const 695 694 { 696 double * const cell_size = World::get()->cell_size;697 695 bool result = true; 698 696 int j =-1; … … 770 768 void molecule::OutputListOfBonds() const 771 769 { 772 DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);770 Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl; 773 771 ActOnAllAtoms (&atom::OutputBondOfAtom ); 774 DoLog(0) && (Log() << Verbose(0) << endl);772 Log() << Verbose(0) << endl; 775 773 }; 776 774 … … 829 827 } 830 828 if ((AtomCount == 0) || (i != AtomCount)) { 831 DoLog(3) && (Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl);829 Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl; 832 830 AtomCount = i; 833 831 … … 845 843 Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name"); 846 844 sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1); 847 DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl);845 Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl; 848 846 i++; 849 847 } 850 848 } else 851 DoLog(3) && (Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl);849 Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl; 852 850 } 853 851 }; … … 909 907 bool result = true; // status of comparison 910 908 911 DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);909 Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl; 912 910 /// first count both their atoms and elements and update lists thereby ... 913 911 //Log() << Verbose(0) << "Counting atoms, updating list" << endl; … … 921 919 if (result) { 922 920 if (AtomCount != OtherMolecule->AtomCount) { 923 DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl);921 Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl; 924 922 result = false; 925 923 } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl; … … 928 926 if (result) { 929 927 if (ElementCount != OtherMolecule->ElementCount) { 930 DoLog(4) && (Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl);928 Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl; 931 929 result = false; 932 930 } else Log() << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl; … … 940 938 } 941 939 if (flag < MAX_ELEMENTS) { 942 DoLog(4) && (Log() << Verbose(4) << "ElementsInMolecule don't match." << endl);940 Log() << Verbose(4) << "ElementsInMolecule don't match." << endl; 943 941 result = false; 944 942 } else Log() << Verbose(4) << "ElementsInMolecule match." << endl; … … 946 944 /// then determine and compare center of gravity for each molecule ... 947 945 if (result) { 948 DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);946 Log() << Verbose(5) << "Calculating Centers of Gravity" << endl; 949 947 DeterminePeriodicCenter(CenterOfGravity); 950 948 OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity); 951 DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: ");949 Log() << Verbose(5) << "Center of Gravity: "; 952 950 CenterOfGravity.Output(); 953 DoLog(0) && (Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ");951 Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: "; 954 952 OtherCenterOfGravity.Output(); 955 DoLog(0) && (Log() << Verbose(0) << endl);953 Log() << Verbose(0) << endl; 956 954 if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) { 957 DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);955 Log() << Verbose(4) << "Centers of gravity don't match." << endl; 958 956 result = false; 959 957 } … … 962 960 /// ... then make a list with the euclidian distance to this center for each atom of both molecules 963 961 if (result) { 964 DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);962 Log() << Verbose(5) << "Calculating distances" << endl; 965 963 Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances"); 966 964 OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances"); … … 969 967 970 968 /// ... sort each list (using heapsort (o(N log N)) from GSL) 971 DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);969 Log() << Verbose(5) << "Sorting distances" << endl; 972 970 PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap"); 973 971 OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap"); … … 975 973 gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles); 976 974 PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap"); 977 DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);975 Log() << Verbose(5) << "Combining Permutation Maps" << endl; 978 976 for(int i=AtomCount;i--;) 979 977 PermutationMap[PermMap[i]] = (int) OtherPermMap[i]; 980 978 981 979 /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all 982 DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);980 Log() << Verbose(4) << "Comparing distances" << endl; 983 981 flag = 0; 984 982 for (int i=0;i<AtomCount;i++) { 985 DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl);983 Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl; 986 984 if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold) 987 985 flag = 1; … … 999 997 } 1000 998 /// return pointer to map if all distances were below \a threshold 1001 DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);999 Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl; 1002 1000 if (result) { 1003 DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);1001 Log() << Verbose(3) << "Result: Equal." << endl; 1004 1002 return PermutationMap; 1005 1003 } else { 1006 DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);1004 Log() << Verbose(3) << "Result: Not equal." << endl; 1007 1005 return NULL; 1008 1006 } … … 1019 1017 { 1020 1018 atom *Walker = NULL, *OtherWalker = NULL; 1021 DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl);1019 Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl; 1022 1020 int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap"); 1023 1021 for (int i=AtomCount;i--;) … … 1026 1024 for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence 1027 1025 AtomicMap[i] = i; 1028 DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl);1026 Log() << Verbose(4) << "Map is trivial." << endl; 1029 1027 } else { 1030 DoLog(4) && (Log() << Verbose(4) << "Map is ");1028 Log() << Verbose(4) << "Map is "; 1031 1029 Walker = start; 1032 1030 while (Walker->next != end) { … … 1045 1043 } 1046 1044 } 1047 DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t");1048 } 1049 DoLog(0) && (Log() << Verbose(0) << endl);1050 } 1051 DoLog(3) && (Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl);1045 Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t"; 1046 } 1047 Log() << Verbose(0) << endl; 1048 } 1049 Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl; 1052 1050 return AtomicMap; 1053 1051 };
Note:
See TracChangeset
for help on using the changeset viewer.