Changeset 6250e5 for molecuilder/src/analysis_bonds.cpp

Timestamp:

Apr 1, 2010, 6:52:09 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Children:

6dd8d3

Parents:

794482

Message:

Enhanced CountHydrogenBridgeBonds() a bit.

• new functions CheckHydrogenBridgeBondAngle() and CalculateAngle()
• CountHydrogenBridgeBonds() now discards when two oxygen show each both a hydrogen (two test cases in CountBondsUnitTest.cpp for this)
• this however also no effect over the criteria by Marry, Rotenberg and Turq.

File:

• molecuilder/src/analysis_bonds.cpp (modified) (3 diffs)

molecuilder/src/analysis_bonds.cpp

@@ -82 +82 @@
 };
 
+/** Calculate the angle between \a *first and \a *origin and \a *second and \a *origin.
+ * \param *first first Vector
+ * \param *origin origin of angle taking
+ * \param *second second Vector
+ * \return angle between \a *first and \a *second, both relative to origin at \a *origin.
+ */
+double CalculateAngle(Vector *first, Vector *central, Vector *second)
+{
+  Vector OHBond;
+  Vector OOBond;
+
+  OHBond.CopyVector(first);
+  OHBond.SubtractVector(central);
+  OOBond.CopyVector(second);
+  OOBond.SubtractVector(central);
+  const double angle = OHBond.Angle(&OOBond);
+  return angle;
+};
+
+/** Checks whether the angle between \a *Oxygen and \a *Hydrogen and \a *Oxygen and \a *OtherOxygen is less than 30 degrees.
+ * Note that distance criterion is not checked.
+ * \param *Oxygen first oxygen atom, bonded to \a *Hydrogen
+ * \param *Hydrogen hydrogen bonded to \a *Oxygen
+ * \param *OtherOxygen other oxygen atom
+ * \return true - angle criteria fulfilled, false - criteria not fulfilled, angle greater than 30 degrees.
+ */
+bool CheckHydrogenBridgeBondAngle(atom *Oxygen, atom *Hydrogen, atom *OtherOxygen)
+{
+  Info FunctionInfo(__func__);
+
+  // check angle
+  if (CalculateAngle(&Hydrogen->x, &Oxygen->x, &OtherOxygen->x) < M_PI*(30./180.)) {
+    return true;
+  } else {
+    return false;
+  }
+};
 
 /** Counts the number of hydrogen bridge bonds.
@@ -91 +128 @@
 int CountHydrogenBridgeBonds(MoleculeListClass *molecules, element * InterfaceElement = NULL)
 {
-  Info FunctionInfo(__func__);
   atom *Walker = NULL;
   atom *Runner = NULL;
-  Vector OHBond;
-  Vector OOBond;
   int count = 0;
+  int OtherHydrogens = 0;
+  double Otherangle = 0.;
   bool InterfaceFlag = false;
+  bool OtherHydrogenFlag = true;
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin();MolWalker != molecules->ListOfMolecules.end(); MolWalker++) {
     Walker = (*MolWalker)->start;
     while (Walker->next != (*MolWalker)->end) {
       Walker = Walker->next;
-      for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != MolWalker; MolRunner++) {
+      for (MoleculeList::const_iterator MolRunner = molecules->ListOfMolecules.begin();MolRunner != molecules->ListOfMolecules.end(); MolRunner++) {
         Runner = (*MolRunner)->start;
         while (Runner->next != (*MolRunner)->end) {
           Runner = Runner->next;
-          if ((Runner != Walker) && (Walker->type->Z  == 8) && (Runner->type->Z  == 8)) {
+          if ((Walker->type->Z  == 8) && (Runner->type->Z  == 8)) {
             // check distance
             const double distance = Runner->x.DistanceSquared(&Walker->x);
             if ((distance > MYEPSILON) && (distance < HBRIDGEDISTANCE*HBRIDGEDISTANCE)) { // distance >0 means  different atoms
+              // on other atom(Runner) we check for bond to interface element and
+              // check that O-O line is not in between the shanks of the two connected hydrogens (Otherangle > 104.5)
+              OtherHydrogenFlag = true;
+              Otherangle = 0.;
+              OtherHydrogens = 0;
               InterfaceFlag = (InterfaceElement == NULL);
-              // on other atom(Runner) we check for bond to interface element
-              for (BondList::const_iterator BondRunner = Runner->ListOfBonds.begin(); ((!InterfaceFlag) && (BondRunner != Runner->ListOfBonds.end())); BondRunner++) {
+              for (BondList::const_iterator BondRunner = Runner->ListOfBonds.begin(); BondRunner != Runner->ListOfBonds.end(); BondRunner++) {
                 atom * const OtherAtom = (*BondRunner)->GetOtherAtom(Runner);
+                // if hydrogen, check angle to be greater(!) than 30 degrees
+                if (OtherAtom->type->Z == 1) {
+                  const double angle = CalculateAngle(&OtherAtom->x, &Runner->x, &Walker->x);
+                  OtherHydrogenFlag = OtherHydrogenFlag && (angle > M_PI*(30./180.) + MYEPSILON);
+                  Otherangle += angle;
+                  OtherHydrogens++;
+                }
                 InterfaceFlag = InterfaceFlag || (OtherAtom->type == InterfaceElement);
               }
-              if (InterfaceFlag) {
+              DoLog(1) && (Log() << Verbose(1) << "Otherangle is " << Otherangle << " for " << OtherHydrogens << " hydrogens." << endl);
+              switch (OtherHydrogens) {
+                case 0:
+                case 1:
+                  break;
+                case 2:
+                  OtherHydrogenFlag = OtherHydrogenFlag && (Otherangle > M_PI*(104.5/180.) + MYEPSILON);
+                  break;
+                default: // 3 or more hydrogens ...
+                  OtherHydrogenFlag = false;
+                  break;
+              }
+              if (InterfaceFlag && OtherHydrogenFlag) {
                 // 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++) {
@@ -122 +182 @@
                   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);
+                    if (CheckHydrogenBridgeBondAngle(Walker, OtherAtom, Runner)) {
+                      DoLog(1) && (Log() << Verbose(1) << Walker->Name << ", " << OtherAtom->Name << " and " << Runner->Name << " has a hydrogen bridge bond with distance " << sqrt(distance) << " and angle " << CalculateAngle(&OtherAtom->x, &Walker->x, &Runner->x)*(180./M_PI) << "." << endl);
                       count++;
                       break;