Changeset 5f612ee for src/molecule.cpp
- Timestamp:
- Apr 27, 2010, 2:25:42 PM (15 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, 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:
- 632bc3
- Parents:
- 13d5a9 (diff), c695c9 (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
Legend:
- Unmodified
- Added
- Removed
-
src/molecule.cpp
r13d5a9 r5f612ee 25 25 #include "tesselation.hpp" 26 26 #include "vector.hpp" 27 #include "World.hpp" 27 28 28 29 /************************************* Functions for class molecule *********************************/ … … 50 51 for(int i=MAX_ELEMENTS;i--;) 51 52 ElementsInMolecule[i] = 0; 52 cell_size[0] = cell_size[2] = cell_size[5]= 20.; 53 cell_size[1] = cell_size[3] = cell_size[4]= 0.; 54 strcpy(name,"none"); 53 strcpy(name,World::getInstance().getDefaultName()); 55 54 }; 56 55 … … 215 214 double *matrix = NULL; 216 215 bond *Binder = NULL; 216 double * const cell_size = World::getInstance().getDomain(); 217 217 218 218 // Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl; … … 250 250 BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1]; 251 251 if (BondRescale == -1) { 252 eLog() << Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;252 DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl); 253 253 return false; 254 254 BondRescale = bondlength; … … 293 293 SecondOtherAtom = (*Runner)->GetOtherAtom(TopOrigin); 294 294 } else { 295 eLog() << Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name;295 DoeLog(2) && (eLog()<< Verbose(2) << "Detected more than four bonds for atom " << TopOrigin->Name); 296 296 } 297 297 } … … 330 330 bondangle = TopOrigin->type->HBondAngle[1]; 331 331 if (bondangle == -1) { 332 eLog() << Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;332 DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl); 333 333 return false; 334 334 bondangle = 0; … … 452 452 break; 453 453 default: 454 eLog() << Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl;454 DoeLog(1) && (eLog()<< Verbose(1) << "BondDegree does not state single, double or triple bond!" << endl); 455 455 AllWentWell = false; 456 456 break; … … 487 487 input = new istringstream(line); 488 488 *input >> NumberOfAtoms; 489 Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;489 DoLog(0) && (Log() << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl); 490 490 getline(xyzfile,line,'\n'); // Read comment 491 Log() << Verbose(1) << "Comment: " << line << endl;491 DoLog(1) && (Log() << Verbose(1) << "Comment: " << line << endl); 492 492 493 493 if (MDSteps == 0) // no atoms yet present … … 505 505 Walker->type = elemente->FindElement(shorthand); 506 506 if (Walker->type == NULL) { 507 eLog() << Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.";507 DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H."); 508 508 Walker->type = elemente->FindElement(1); 509 509 } … … 532 532 molecule *molecule::CopyMolecule() 533 533 { 534 molecule *copy = new molecule(elemente);534 molecule *copy = World::getInstance().createMolecule(); 535 535 atom *LeftAtom = NULL, *RightAtom = NULL; 536 536 … … 575 575 */ 576 576 molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const { 577 molecule *copy = new molecule(elemente);577 molecule *copy = World::getInstance().createMolecule(); 578 578 579 579 ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped ); … … 601 601 add(Binder, last); 602 602 } else { 603 eLog() << Verbose(1) << "Could not add bond between " << atom1->Name << " and " << atom2->Name << " as one or both are not present in the molecule." << endl;603 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); 604 604 } 605 605 return Binder; … … 613 613 bool molecule::RemoveBond(bond *pointer) 614 614 { 615 // eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;615 //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl); 616 616 pointer->leftatom->RegisterBond(pointer); 617 617 pointer->rightatom->RegisterBond(pointer); … … 627 627 bool molecule::RemoveBonds(atom *BondPartner) 628 628 { 629 // eLog() << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;629 //DoeLog(1) && (eLog()<< Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl); 630 630 BondList::const_iterator ForeRunner; 631 631 while (!BondPartner->ListOfBonds.empty()) { … … 661 661 void molecule::SetBoxDimension(Vector *dim) 662 662 { 663 double * const cell_size = World::getInstance().getDomain(); 663 664 cell_size[0] = dim->x[0]; 664 665 cell_size[1] = 0.; … … 679 680 AtomCount--; 680 681 } else 681 eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;682 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); 682 683 if (ElementsInMolecule[pointer->type->Z] == 0) // was last atom of this element? 683 684 ElementCount--; … … 697 698 ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element 698 699 else 699 eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;700 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); 700 701 if (ElementsInMolecule[pointer->type->Z] == 0) // was last atom of this element? 701 702 ElementCount--; … … 722 723 return walker; 723 724 } else { 724 Log() << Verbose(0) << "Atom not found in list." << endl;725 DoLog(0) && (Log() << Verbose(0) << "Atom not found in list." << endl); 725 726 return NULL; 726 727 } … … 738 739 //mol->Output((ofstream *)&cout); 739 740 //Log() << Verbose(0) << "===============================================" << endl; 740 Log() << Verbose(0) << text;741 DoLog(0) && (Log() << Verbose(0) << text); 741 742 cin >> No; 742 743 ion = this->FindAtom(No); … … 751 752 bool molecule::CheckBounds(const Vector *x) const 752 753 { 754 double * const cell_size = World::getInstance().getDomain(); 753 755 bool result = true; 754 756 int j =-1; … … 826 828 void molecule::OutputListOfBonds() const 827 829 { 828 Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl;830 DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl); 829 831 ActOnAllAtoms (&atom::OutputBondOfAtom ); 830 Log() << Verbose(0) << endl;832 DoLog(0) && (Log() << Verbose(0) << endl); 831 833 }; 832 834 … … 885 887 } 886 888 if ((AtomCount == 0) || (i != AtomCount)) { 887 Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;889 DoLog(3) && (Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl); 888 890 AtomCount = i; 889 891 … … 901 903 Walker->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name"); 902 904 sprintf(Walker->Name, "%2s%02d", Walker->type->symbol, Walker->nr+1); 903 Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl;905 DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << Walker->nr << " " << Walker->Name << "." << endl); 904 906 i++; 905 907 } 906 908 } else 907 Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;909 DoLog(3) && (Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl); 908 910 } 909 911 }; … … 965 967 bool result = true; // status of comparison 966 968 967 Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;969 DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl); 968 970 /// first count both their atoms and elements and update lists thereby ... 969 971 //Log() << Verbose(0) << "Counting atoms, updating list" << endl; … … 977 979 if (result) { 978 980 if (AtomCount != OtherMolecule->AtomCount) { 979 Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;981 DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl); 980 982 result = false; 981 983 } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl; … … 984 986 if (result) { 985 987 if (ElementCount != OtherMolecule->ElementCount) { 986 Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl;988 DoLog(4) && (Log() << Verbose(4) << "ElementCount don't match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl); 987 989 result = false; 988 990 } else Log() << Verbose(4) << "ElementCount match: " << ElementCount << " == " << OtherMolecule->ElementCount << endl; … … 996 998 } 997 999 if (flag < MAX_ELEMENTS) { 998 Log() << Verbose(4) << "ElementsInMolecule don't match." << endl;1000 DoLog(4) && (Log() << Verbose(4) << "ElementsInMolecule don't match." << endl); 999 1001 result = false; 1000 1002 } else Log() << Verbose(4) << "ElementsInMolecule match." << endl; … … 1002 1004 /// then determine and compare center of gravity for each molecule ... 1003 1005 if (result) { 1004 Log() << Verbose(5) << "Calculating Centers of Gravity" << endl;1006 DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl); 1005 1007 DeterminePeriodicCenter(CenterOfGravity); 1006 1008 OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity); 1007 Log() << Verbose(5) << "Center of Gravity: ";1009 DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: "); 1008 1010 CenterOfGravity.Output(); 1009 Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: ";1011 DoLog(0) && (Log() << Verbose(0) << endl << Verbose(5) << "Other Center of Gravity: "); 1010 1012 OtherCenterOfGravity.Output(); 1011 Log() << Verbose(0) << endl;1013 DoLog(0) && (Log() << Verbose(0) << endl); 1012 1014 if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) { 1013 Log() << Verbose(4) << "Centers of gravity don't match." << endl;1015 DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl); 1014 1016 result = false; 1015 1017 } … … 1018 1020 /// ... then make a list with the euclidian distance to this center for each atom of both molecules 1019 1021 if (result) { 1020 Log() << Verbose(5) << "Calculating distances" << endl;1022 DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl); 1021 1023 Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances"); 1022 1024 OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances"); … … 1025 1027 1026 1028 /// ... sort each list (using heapsort (o(N log N)) from GSL) 1027 Log() << Verbose(5) << "Sorting distances" << endl;1029 DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl); 1028 1030 PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap"); 1029 1031 OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap"); … … 1031 1033 gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles); 1032 1034 PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap"); 1033 Log() << Verbose(5) << "Combining Permutation Maps" << endl;1035 DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl); 1034 1036 for(int i=AtomCount;i--;) 1035 1037 PermutationMap[PermMap[i]] = (int) OtherPermMap[i]; 1036 1038 1037 1039 /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all 1038 Log() << Verbose(4) << "Comparing distances" << endl;1040 DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl); 1039 1041 flag = 0; 1040 1042 for (int i=0;i<AtomCount;i++) { 1041 Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl;1043 DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl); 1042 1044 if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold) 1043 1045 flag = 1; … … 1055 1057 } 1056 1058 /// return pointer to map if all distances were below \a threshold 1057 Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;1059 DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl); 1058 1060 if (result) { 1059 Log() << Verbose(3) << "Result: Equal." << endl;1061 DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl); 1060 1062 return PermutationMap; 1061 1063 } else { 1062 Log() << Verbose(3) << "Result: Not equal." << endl;1064 DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl); 1063 1065 return NULL; 1064 1066 } … … 1075 1077 { 1076 1078 atom *Walker = NULL, *OtherWalker = NULL; 1077 Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;1079 DoLog(3) && (Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl); 1078 1080 int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap"); 1079 1081 for (int i=AtomCount;i--;) … … 1082 1084 for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence 1083 1085 AtomicMap[i] = i; 1084 Log() << Verbose(4) << "Map is trivial." << endl;1086 DoLog(4) && (Log() << Verbose(4) << "Map is trivial." << endl); 1085 1087 } else { 1086 Log() << Verbose(4) << "Map is ";1088 DoLog(4) && (Log() << Verbose(4) << "Map is "); 1087 1089 Walker = start; 1088 1090 while (Walker->next != end) { … … 1101 1103 } 1102 1104 } 1103 Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t";1104 } 1105 Log() << Verbose(0) << endl;1106 } 1107 Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl;1105 DoLog(0) && (Log() << Verbose(0) << AtomicMap[Walker->nr] << "\t"); 1106 } 1107 DoLog(0) && (Log() << Verbose(0) << endl); 1108 } 1109 DoLog(3) && (Log() << Verbose(3) << "End of GetFatherAtomicMap." << endl); 1108 1110 return AtomicMap; 1109 1111 };
Note:
See TracChangeset
for help on using the changeset viewer.