Changeset 822f01
- Timestamp:
- Aug 6, 2010, 2:09:12 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:
- e588312
- Parents:
- e5c0a1
- Location:
- src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
src/molecule.cpp
re5c0a1 r822f01 983 983 }; 984 984 985 /** Determines whether two molecules actually contain the same atoms and coordination.986 * \param *out output stream for debugging987 * \param *OtherMolecule the molecule to compare this one to988 * \param threshold upper limit of difference when comparing the coordination.989 * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)990 */991 int * molecule::IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold)992 {993 int flag;994 double *Distances = NULL, *OtherDistances = NULL;995 Vector CenterOfGravity, OtherCenterOfGravity;996 size_t *PermMap = NULL, *OtherPermMap = NULL;997 int *PermutationMap = NULL;998 bool result = true; // status of comparison999 1000 DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);1001 /// first count both their atoms and elements and update lists thereby ...1002 //Log() << Verbose(0) << "Counting atoms, updating list" << endl;1003 1004 /// ... and compare:1005 /// -# AtomCount1006 if (result) {1007 if (getAtomCount() != OtherMolecule->getAtomCount()) {1008 DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl);1009 result = false;1010 } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;1011 }1012 /// -# Formula1013 if (result) {1014 if (formula != OtherMolecule->formula) {1015 DoLog(4) && (Log() << Verbose(4) << "Formulas don't match: " << formula << " == " << OtherMolecule->formula << endl);1016 result = false;1017 } else Log() << Verbose(4) << "Formulas match: " << formula << " == " << OtherMolecule->formula << endl;1018 }1019 /// then determine and compare center of gravity for each molecule ...1020 if (result) {1021 DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);1022 DeterminePeriodicCenter(CenterOfGravity);1023 OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);1024 DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: " << CenterOfGravity << endl);1025 DoLog(5) && (Log() << Verbose(5) << "Other Center of Gravity: " << OtherCenterOfGravity << endl);1026 if (CenterOfGravity.DistanceSquared(OtherCenterOfGravity) > threshold*threshold) {1027 DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);1028 result = false;1029 }1030 }1031 1032 /// ... then make a list with the euclidian distance to this center for each atom of both molecules1033 if (result) {1034 DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);1035 Distances = new double[getAtomCount()];1036 OtherDistances = new double[getAtomCount()];1037 SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);1038 SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);1039 for(int i=0;i<getAtomCount();i++) {1040 Distances[i] = 0.;1041 OtherDistances[i] = 0.;1042 }1043 1044 /// ... sort each list (using heapsort (o(N log N)) from GSL)1045 DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);1046 PermMap = new size_t[getAtomCount()];1047 OtherPermMap = new size_t[getAtomCount()];1048 for(int i=0;i<getAtomCount();i++) {1049 PermMap[i] = 0;1050 OtherPermMap[i] = 0;1051 }1052 gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);1053 gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);1054 PermutationMap = new int[getAtomCount()];1055 for(int i=0;i<getAtomCount();i++)1056 PermutationMap[i] = 0;1057 DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);1058 for(int i=getAtomCount();i--;)1059 PermutationMap[PermMap[i]] = (int) OtherPermMap[i];1060 1061 /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all1062 DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);1063 flag = 0;1064 for (int i=0;i<getAtomCount();i++) {1065 DoLog(5) && (Log() << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " << threshold << endl);1066 if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)1067 flag = 1;1068 }1069 1070 // free memory1071 delete[](PermMap);1072 delete[](OtherPermMap);1073 delete[](Distances);1074 delete[](OtherDistances);1075 if (flag) { // if not equal1076 delete[](PermutationMap);1077 result = false;1078 }1079 }1080 /// return pointer to map if all distances were below \a threshold1081 DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);1082 if (result) {1083 DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);1084 return PermutationMap;1085 } else {1086 DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);1087 return NULL;1088 }1089 };1090 1091 985 /** Returns an index map for two father-son-molecules. 1092 986 * The map tells which atom in this molecule corresponds to which one in the other molecul with their fathers. -
src/molecule.hpp
re5c0a1 r822f01 337 337 338 338 // Recognize doubly appearing molecules in a list of them 339 int * IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold);340 339 int * GetFatherSonAtomicMap(molecule *OtherMolecule); 341 340
Note:
See TracChangeset
for help on using the changeset viewer.