Changeset b5c53d for src/molecule.cpp
- Timestamp:
- Aug 9, 2010, 6:33:29 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, 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:
- 353e82
- Parents:
- 458447 (diff), e588312 (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) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/molecule.cpp
r458447 rb5c53d 218 218 if(pointer->getName() == "Unknown"){ 219 219 stringstream sstr; 220 sstr << pointer->getType()-> symbol<< pointer->nr+1;220 sstr << pointer->getType()->getSymbol() << pointer->nr+1; 221 221 pointer->setName(sstr.str()); 222 222 } … … 972 972 NoNonHydrogen++; 973 973 stringstream sstr; 974 sstr << (*iter)->getType()-> symbol<< (*iter)->nr+1;974 sstr << (*iter)->getType()->getSymbol() << (*iter)->nr+1; 975 975 (*iter)->setName(sstr.str()); 976 976 DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl); … … 978 978 } 979 979 return res; 980 };981 982 /** Determines whether two molecules actually contain the same atoms and coordination.983 * \param *out output stream for debugging984 * \param *OtherMolecule the molecule to compare this one to985 * \param threshold upper limit of difference when comparing the coordination.986 * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)987 */988 int * molecule::IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold)989 {990 int flag;991 double *Distances = NULL, *OtherDistances = NULL;992 Vector CenterOfGravity, OtherCenterOfGravity;993 size_t *PermMap = NULL, *OtherPermMap = NULL;994 int *PermutationMap = NULL;995 bool result = true; // status of comparison996 997 DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);998 /// first count both their atoms and elements and update lists thereby ...999 //Log() << Verbose(0) << "Counting atoms, updating list" << endl;1000 1001 /// ... and compare:1002 /// -# AtomCount1003 if (result) {1004 if (getAtomCount() != OtherMolecule->getAtomCount()) {1005 DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl);1006 result = false;1007 } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;1008 }1009 /// -# Formula1010 if (result) {1011 if (formula != OtherMolecule->formula) {1012 DoLog(4) && (Log() << Verbose(4) << "Formulas don't match: " << formula << " == " << OtherMolecule->formula << endl);1013 result = false;1014 } else Log() << Verbose(4) << "Formulas match: " << formula << " == " << OtherMolecule->formula << endl;1015 }1016 /// then determine and compare center of gravity for each molecule ...1017 if (result) {1018 DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);1019 DeterminePeriodicCenter(CenterOfGravity);1020 OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);1021 DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: " << CenterOfGravity << endl);1022 DoLog(5) && (Log() << Verbose(5) << "Other Center of Gravity: " << OtherCenterOfGravity << endl);1023 if (CenterOfGravity.DistanceSquared(OtherCenterOfGravity) > threshold*threshold) {1024 DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);1025 result = false;1026 }1027 }1028 1029 /// ... then make a list with the euclidian distance to this center for each atom of both molecules1030 if (result) {1031 DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);1032 Distances = new double[getAtomCount()];1033 OtherDistances = new double[getAtomCount()];1034 SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);1035 SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);1036 for(int i=0;i<getAtomCount();i++) {1037 Distances[i] = 0.;1038 OtherDistances[i] = 0.;1039 }1040 1041 /// ... sort each list (using heapsort (o(N log N)) from GSL)1042 DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);1043 PermMap = new size_t[getAtomCount()];1044 OtherPermMap = new size_t[getAtomCount()];1045 for(int i=0;i<getAtomCount();i++) {1046 PermMap[i] = 0;1047 OtherPermMap[i] = 0;1048 }1049 gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);1050 gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);1051 PermutationMap = new int[getAtomCount()];1052 for(int i=0;i<getAtomCount();i++)1053 PermutationMap[i] = 0;1054 DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);1055 for(int i=getAtomCount();i--;)1056 PermutationMap[PermMap[i]] = (int) OtherPermMap[i];1057 1058 /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all1059 DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);1060 flag = 0;1061 for (int i=0;i<getAtomCount();i++) {1062 DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl);1063 if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)1064 flag = 1;1065 }1066 1067 // free memory1068 delete[](PermMap);1069 delete[](OtherPermMap);1070 delete[](Distances);1071 delete[](OtherDistances);1072 if (flag) { // if not equal1073 delete[](PermutationMap);1074 result = false;1075 }1076 }1077 /// return pointer to map if all distances were below \a threshold1078 DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);1079 if (result) {1080 DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);1081 return PermutationMap;1082 } else {1083 DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);1084 return NULL;1085 }1086 980 }; 1087 981
Note:
See TracChangeset
for help on using the changeset viewer.
