Changeset dfe8ef for src/bondgraph.cpp

Timestamp:

Mar 30, 2010, 6:52:14 PM (15 years ago)

Author:

Frederik Heber <heber@…>

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:

62c9c0

Parents:

f18185

git-author:

Frederik Heber <heber@…> (03/30/10 18:49:01)

git-committer:

Frederik Heber <heber@…> (03/30/10 18:52:14)

Message:

BUGFIX: CountHydrogenBridgeBonds() was half-broken, CountBondsOfThree() double-counted.

• CountBondsOfThree() was counting O-H molecules as having a H-O-H and also a O-H-O triple bond because of a missing break.
• CountHydrogenBridgeBonds() was broken in several places:
  • Checking for Hydrogen on the other oxygen atom is not part of the criteria by Marry, Rotenberg and Turq.
  • HydrogenFlag could have been set by previous atom making next hydrogen also a h-bridge bond.

Found thanks to newly implemented unit test CountBondsUnitTest.

Signed-off-by: Frederik Heber <heber@…>

File:

• src/bondgraph.cpp (modified) (3 diffs)

src/bondgraph.cpp

@@ -183 +183 @@
   atom *Walker = NULL;
   atom *Runner = NULL;
-  atom *Hydrogen = NULL;
-  atom *OtherHydrogen = NULL;
   Vector OHBond;
   Vector OOBond;
   int count = 0;
-  bool HydrogenFlag = false;
-  bool OtherHydrogenFlag = false;
   bool InterfaceFlag = false;
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
@@ -202 +198 @@
             // check distance
             const double distance = Runner->x.DistanceSquared(&Walker->x);
-            if (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE) {
-              // get hydrogen, check for InterfaceElement
-              HydrogenFlag = false;
-              OtherHydrogenFlag = false;
+            if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means  different atoms
               InterfaceFlag = (InterfaceElement == NULL);
               // on other atom(Runner) we check for bond to interface element
-              for (BondList::const_iterator BondRunner = Runner->ListOfBonds.begin(); BondRunner != Runner->ListOfBonds.end(); BondRunner++) {
+              for (BondList::const_iterator BondRunner = Runner->ListOfBonds.begin(); ((!InterfaceFlag) && (BondRunner != Runner->ListOfBonds.end())); BondRunner++) {
                 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Runner);
-                if (!OtherHydrogenFlag && (OtherAtom->type->Z == 1)) {
-                  OtherHydrogen = OtherAtom;
-                  OtherHydrogenFlag = true;
-                }
                 InterfaceFlag = InterfaceFlag || (OtherAtom->type == InterfaceElement);
               }
-              // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
-              for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {
-                atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
-                if (!HydrogenFlag && (OtherAtom->type->Z == 1)) {
-                  Hydrogen = OtherAtom;
-                  HydrogenFlag = true;
-                }
-              }
-              if (InterfaceFlag && HydrogenFlag) {
-                if ((Walker->nr < Runner->nr) || (!OtherHydrogenFlag)) {
-                  // check angle
-                  OHBond.CopyVector(&Hydrogen->x);
-                  OHBond.SubtractVector(&Walker->x);
-                  OOBond.CopyVector(&Runner->x);
-                  OOBond.SubtractVector(&Walker->x);
-                  const double angle = OHBond.Angle(&OOBond);
-                  if (angle < M_PI*(30./180.)) {
-                    DoLog(1) && (Log() << Verbose(1) << Walker->Name << ", " << Hydrogen->Name << " and " << Runner->Name << " have a hydrogen bridge bond with " << sqrt(distance) << " and at angle " << (180./M_PI)*angle << " degrees." << endl);
-                    count++;
+              if (InterfaceFlag) {
+                // on this element (Walker) we check for bond to hydrogen, i.e. part of water molecule
+                for (BondList::const_iterator BondRunner = Walker->ListOfBonds.begin(); BondRunner != Walker->ListOfBonds.end(); BondRunner++) {
+                  atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
+                  if (OtherAtom->type->Z == 1) {
+                    // check angle
+                    OHBond.CopyVector(&OtherAtom->x);
+                    OHBond.SubtractVector(&Walker->x);
+                    OOBond.CopyVector(&Runner->x);
+                    OOBond.SubtractVector(&Walker->x);
+                    const double angle = OHBond.Angle(&OOBond);
+                    if (angle < M_PI*(30./180.)) {
+                      DoLog(1) && (Log() << Verbose(1) << Walker->Name << ", " << OtherAtom->Name << " and " << Runner->Name << " have a hydrogen bridge bond with " << sqrt(distance) << " and at angle " << (180./M_PI)*angle << " degrees." << endl);
+                      count++;
+                      break;
+                    }
                   }
                 }
@@ -305 +291 @@
           atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Walker);
           for (int i=0;i<2;i++)
-            MatchFlag[i] = MatchFlag[i] || (OtherAtom->type == ElementArray[i]);
+            if ((!MatchFlag[i]) && (OtherAtom->type == ElementArray[i])) {
+              MatchFlag[i] = true;
+              break;  // each bonding atom can match at most one element we are looking for
+            }
         }
         result = true;
-        for (int i=0;i<2;i++)
+        for (int i=0;i<2;i++) // gather results
           result = result && MatchFlag[i];
-        if (result) {
+        if (result) { // check results
           count++;
           DoLog(1) && (Log() << Verbose(1) << first->name << "-" << second->name << "-" << third->name << " bond found at " << *Walker << "." << endl);