Changeset a7b761b for src/molecule.cpp
- Timestamp:
- May 27, 2010, 10:46:54 AM (16 years ago)
- Branches:
- Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
- Children:
- 1024cb
- Parents:
- 8f215d (diff), 05a97c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the(diff)links above to see all the changes relative to each parent. - File:
-
- 1 edited
-
src/molecule.cpp (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/molecule.cpp
r8f215d ra7b761b 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), 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),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), 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 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 42 AtomCount(this,boost::bind(&molecule::doCountAtoms,this)), last_atom(0), InternalPointer(begin()) 43 { 50 44 // init bond chain list 51 45 link(first,last); … … 69 63 delete(first); 70 64 delete(last); 71 end->getWorld()->destroyAtom(end);72 start->getWorld()->destroyAtom(start);73 65 }; 74 66 … … 81 73 const std::string molecule::getName(){ 82 74 return std::string(name); 75 } 76 77 int molecule::getAtomCount() const{ 78 return *AtomCount; 83 79 } 84 80 … … 104 100 stringstream sstr; 105 101 periodentafel *periode = World::getInstance().getPeriode(); 106 for (atom *Walker = start; Walker != end; Walker = Walker->next) {107 counts[ Walker->type->getNumber()]++;102 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 103 counts[(*iter)->type->getNumber()]++; 108 104 } 109 105 std::map<atomicNumber_t,unsigned int>::reverse_iterator iter; … … 115 111 } 116 112 113 /************************** Access to the List of Atoms ****************/ 114 115 116 molecule::iterator molecule::begin(){ 117 return molecule::iterator(atoms.begin(),this); 118 } 119 120 molecule::const_iterator molecule::begin() const{ 121 return atoms.begin(); 122 } 123 124 molecule::iterator molecule::end(){ 125 return molecule::iterator(atoms.end(),this); 126 } 127 128 molecule::const_iterator molecule::end() const{ 129 return atoms.end(); 130 } 131 132 bool molecule::empty() const 133 { 134 return (begin() == end()); 135 } 136 137 size_t molecule::size() const 138 { 139 size_t counter = 0; 140 for (molecule::const_iterator iter = begin(); iter != end (); ++iter) 141 counter++; 142 return counter; 143 } 144 145 molecule::const_iterator molecule::erase( const_iterator loc ) 146 { 147 molecule::const_iterator iter = loc; 148 iter--; 149 atoms.erase( loc ); 150 return iter; 151 } 152 153 molecule::const_iterator molecule::erase( atom *& key ) 154 { 155 cout << "trying to erase atom" << endl; 156 molecule::const_iterator iter = find(key); 157 if (iter != end()){ 158 // remove this position and step forward (post-increment) 159 atoms.erase( iter++ ); 160 } 161 return iter; 162 } 163 164 molecule::const_iterator molecule::find ( atom *& key ) const 165 { 166 return atoms.find( key ); 167 } 168 169 pair<molecule::iterator,bool> molecule::insert ( atom * const key ) 170 { 171 pair<atomSet::iterator,bool> res = atoms.insert(key); 172 return pair<iterator,bool>(iterator(res.first,this),res.second); 173 } 117 174 118 175 /** Adds given atom \a *pointer from molecule list. … … 123 180 bool molecule::AddAtom(atom *pointer) 124 181 { 125 bool retval = false;126 182 OBSERVE; 127 183 if (pointer != NULL) { 128 184 pointer->sort = &pointer->nr; 129 pointer->nr = last_atom++; // increase number within molecule130 AtomCount++;131 185 if (pointer->type != NULL) { 132 186 if (ElementsInMolecule[pointer->type->Z] == 0) … … 141 195 } 142 196 } 143 retval = add(pointer, end);144 } 145 return retval;197 insert(pointer); 198 } 199 return true; 146 200 }; 147 201 … … 157 211 if (pointer != NULL) { 158 212 atom *walker = pointer->clone(); 159 stringstream sstr; 160 sstr << pointer->getName(); 161 walker->setName(sstr.str()); 213 walker->setName(pointer->getName()); 162 214 walker->nr = last_atom++; // increase number within molecule 163 add(walker, end);215 insert(walker); 164 216 if ((pointer->type != NULL) && (pointer->type->Z != 1)) 165 217 NoNonHydrogen++; 166 AtomCount++;167 218 retval=walker; 168 219 } … … 574 625 575 626 // copy values 576 copy->CountAtoms();577 627 copy->CountElements(); 578 628 if (first->next != last) { // if adjaceny list is present … … 609 659 { 610 660 bond *Binder = NULL; 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 } 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 621 675 return Binder; 622 676 }; … … 692 746 bool molecule::RemoveAtom(atom *pointer) 693 747 { 748 ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom()."); 749 OBSERVE; 694 750 if (ElementsInMolecule[pointer->type->Z] != 0) { // this would indicate an error 695 751 ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element 696 AtomCount--;697 752 } else 698 753 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); … … 700 755 ElementCount--; 701 756 RemoveBonds(pointer); 702 return remove(pointer, start, end); 757 erase(pointer); 758 return true; 703 759 }; 704 760 … … 717 773 if (ElementsInMolecule[pointer->type->Z] == 0) // was last atom of this element? 718 774 ElementCount--; 719 unlink(pointer);775 erase(pointer); 720 776 return true; 721 777 }; … … 726 782 bool molecule::CleanupMolecule() 727 783 { 728 return (cleanup(first,last) && cleanup(start,end)); 784 for (molecule::iterator iter = begin(); !empty(); iter = begin()) 785 erase(iter); 786 return (cleanup(first,last)); 729 787 }; 730 788 … … 733 791 * \return pointer to atom or NULL 734 792 */ 735 atom * molecule::FindAtom(int Nr) const{ 736 atom * walker = find(&Nr, start,end); 737 if (walker != NULL) { 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()) { 738 800 //Log() << Verbose(0) << "Found Atom Nr. " << walker->nr << endl; 739 return walker;801 return (*iter); 740 802 } else { 741 803 DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl); … … 867 929 now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time' 868 930 for (int step=0;step<MDSteps;step++) { 869 *output << AtomCount<< "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);931 *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now); 870 932 ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step ); 871 933 } … … 884 946 if (output != NULL) { 885 947 now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time' 886 *output << AtomCount<< "\n\tCreated by molecuilder on " << ctime(&now);948 *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now); 887 949 ActOnAllAtoms( &atom::OutputXYZLine, output ); 888 950 return true; … … 894 956 * \param *out output stream for debugging 895 957 */ 896 void molecule::CountAtoms() 897 { 958 int molecule::doCountAtoms() 959 { 960 int res = size(); 898 961 int i = 0; 899 atom *Walker = start; 900 while (Walker->next != end) { 901 Walker = Walker->next; 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; 902 971 i++; 903 972 } 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 } 973 return res; 927 974 }; 928 975 … … 986 1033 /// first count both their atoms and elements and update lists thereby ... 987 1034 //Log() << Verbose(0) << "Counting atoms, updating list" << endl; 988 CountAtoms();989 OtherMolecule->CountAtoms();990 1035 CountElements(); 991 1036 OtherMolecule->CountElements(); … … 994 1039 /// -# AtomCount 995 1040 if (result) { 996 if ( AtomCount != OtherMolecule->AtomCount) {997 DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount<< endl);1041 if (getAtomCount() != OtherMolecule->getAtomCount()) { 1042 DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl); 998 1043 result = false; 999 } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount<< endl;1044 } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl; 1000 1045 } 1001 1046 /// -# ElementCount … … 1034 1079 if (result) { 1035 1080 DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl); 1036 Distances = Calloc<double>( AtomCount, "molecule::IsEqualToWithinThreshold: Distances");1037 OtherDistances = Calloc<double>( AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");1081 Distances = Calloc<double>(getAtomCount(), "molecule::IsEqualToWithinThreshold: Distances"); 1082 OtherDistances = Calloc<double>(getAtomCount(), "molecule::IsEqualToWithinThreshold: OtherDistances"); 1038 1083 SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity); 1039 1084 SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity); … … 1041 1086 /// ... sort each list (using heapsort (o(N log N)) from GSL) 1042 1087 DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl); 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");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"); 1048 1093 DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl); 1049 for(int i= AtomCount;i--;)1094 for(int i=getAtomCount();i--;) 1050 1095 PermutationMap[PermMap[i]] = (int) OtherPermMap[i]; 1051 1096 … … 1053 1098 DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl); 1054 1099 flag = 0; 1055 for (int i=0;i< AtomCount;i++) {1100 for (int i=0;i<getAtomCount();i++) { 1056 1101 DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl); 1057 1102 if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold) … … 1089 1134 int * molecule::GetFatherSonAtomicMap(molecule *OtherMolecule) 1090 1135 { 1091 atom *Walker = NULL, *OtherWalker = NULL;1092 1136 DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl); 1093 int *AtomicMap = Malloc<int>( AtomCount, "molecule::GetAtomicMap: *AtomicMap");1094 for (int i= AtomCount;i--;)1137 int *AtomicMap = Malloc<int>(getAtomCount(), "molecule::GetAtomicMap: *AtomicMap"); 1138 for (int i=getAtomCount();i--;) 1095 1139 AtomicMap[i] = -1; 1096 1140 if (OtherMolecule == this) { // same molecule 1097 for (int i= AtomCount;i--;) // no need as -1 means already that there is trivial correspondence1141 for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence 1098 1142 AtomicMap[i] = i; 1099 1143 DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl); 1100 1144 } else { 1101 1145 DoLog(4) && (Log() << Verbose(4) << "Map is "); 1102 Walker = start; 1103 while (Walker->next != end) { 1104 Walker = Walker->next; 1105 if (Walker->father == NULL) { 1106 AtomicMap[Walker->nr] = -2; 1146 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 1147 if ((*iter)->father == NULL) { 1148 AtomicMap[(*iter)->nr] = -2; 1107 1149 } else { 1108 OtherWalker = OtherMolecule->start; 1109 while (OtherWalker->next != OtherMolecule->end) { 1110 OtherWalker = OtherWalker->next; 1150 for (molecule::const_iterator runner = OtherMolecule->begin(); runner != OtherMolecule->end(); ++runner) { 1111 1151 //for (int i=0;i<AtomCount;i++) { // search atom 1112 1152 //for (int j=0;j<OtherMolecule->AtomCount;j++) { 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;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; 1116 1156 } 1117 1157 } 1118 DoLog(0) && (Log() << Verbose(0) << AtomicMap[ Walker->nr] << "\t");1158 DoLog(0) && (Log() << Verbose(0) << AtomicMap[(*iter)->nr] << "\t"); 1119 1159 } 1120 1160 DoLog(0) && (Log() << Verbose(0) << endl); … … 1150 1190 void molecule::SetIndexedArrayForEachAtomTo ( atom **array, int ParticleInfo::*index) const 1151 1191 { 1152 atom *Walker = start; 1153 while (Walker->next != end) { 1154 Walker = Walker->next; 1155 array[(Walker->*index)] = Walker; 1192 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { 1193 array[((*iter)->*index)] = (*iter); 1156 1194 } 1157 1195 };
Note:
See TracChangeset
for help on using the changeset viewer.
