--- trunk/src/primitives/Molecule.cpp	2013/07/19 21:25:45	1908
+++ trunk/src/primitives/Molecule.cpp	2015/03/07 21:41:51	2071
@@ -53,13 +53,15 @@
 #include "primitives/Molecule.hpp"
 #include "utils/MemoryUtils.hpp"
 #include "utils/simError.h"
+#include "utils/StringUtils.hpp"
 
 namespace OpenMD {
   Molecule::Molecule(int stampId, int globalIndex, const std::string& molName, 
-                     int region) : stampId_(stampId),
+                     int region) : 
                                    globalIndex_(globalIndex), 
-                                   moleculeName_(molName), 
-                                   region_(region),
+                                   stampId_(stampId),
+				   region_(region),
+				   moleculeName_(molName),
                                    constrainTotalCharge_(false) {
   }
   
@@ -145,9 +147,13 @@ namespace OpenMD {
     
     std::set<Atom*> rigidAtoms;
     Atom* atom;
-    AtomIterator ai;
+    Atom* atom1;
+    Atom* atom2;
+    AtomIterator ai, aj;
     RigidBody* rb;
     RigidBodyIterator rbIter;
+    Bond* bond;
+    BondIterator bi;
 
     // Get list of all the atoms that are part of rigid bodies
 
@@ -182,8 +188,89 @@ namespace OpenMD {
         fluctuatingCharges_.push_back( atom );      
     }
 
-  }
 
+    // find the electronegative atoms and add them to the
+    // hBondAcceptors_ vector:
+    
+    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
+      AtomType* at = atom->getAtomType();
+      // get the chain of base types for this atom type:
+      std::vector<AtomType*> ayb = at->allYourBase();
+      // use the last type in the chain of base types for the name:
+      std::string bn = UpperCase(ayb[ayb.size()-1]->getName());
+
+        if (bn.compare("O")==0 || bn.compare("N")==0
+            || bn.compare("F")==0) 
+          hBondAcceptors_.push_back( atom );
+      
+    }
+    
+    // find electronegative atoms that are either bonded to
+    // hydrogens or are present in the same rigid bodies:
+    
+    for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) {
+      Atom* atom1 = bond->getAtomA();
+      Atom* atom2 = bond->getAtomB();
+      AtomType* at1 = atom1->getAtomType();
+      AtomType* at2 = atom1->getAtomType();
+      // get the chain of base types for this atom type:
+      std::vector<AtomType*> ayb1 = at1->allYourBase();
+      std::vector<AtomType*> ayb2 = at2->allYourBase();
+      // use the last type in the chain of base types for the name:
+      std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
+      std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
+      
+      if (bn1.compare("H")==0) {
+        if (bn2.compare("O")==0 || bn2.compare("N")==0
+            || bn2.compare("F")==0) {
+          HBondDonor* donor = new HBondDonor();
+          donor->donorAtom = atom2;
+          donor->donatedHydrogen = atom1;
+          hBondDonors_.push_back( donor );
+        }
+      }
+      if (bn2.compare("H")==0) {
+        if (bn1.compare("O")==0 || bn1.compare("N")==0
+            || bn1.compare("F")==0) {
+          HBondDonor* donor = new HBondDonor();
+          donor->donorAtom = atom1;
+          donor->donatedHydrogen = atom2;
+            hBondDonors_.push_back( donor );
+        }
+      } 
+    }
+    
+    for (rb = beginRigidBody(rbIter); rb != NULL;
+         rb = nextRigidBody(rbIter)) {
+      for(atom1 = rb->beginAtom(ai); atom1 != NULL;
+          atom1 = rb->nextAtom(ai)) {
+        AtomType* at1 = atom1->getAtomType();
+        // get the chain of base types for this atom type:
+        std::vector<AtomType*> ayb1 = at1->allYourBase();
+        // use the last type in the chain of base types for the name:
+        std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
+        
+        if (bn1.compare("O")==0 || bn1.compare("N")==0
+            || bn1.compare("F")==0) {
+          for(atom2 = rb->beginAtom(aj); atom2 != NULL;
+              atom2 = rb->nextAtom(aj)) {
+            AtomType* at2 = atom2->getAtomType();
+            // get the chain of base types for this atom type:              
+            std::vector<AtomType*> ayb2 = at2->allYourBase();
+            // use the last type in the chain of base types for the name:
+            std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
+            if (bn2.compare("H")==0) {              
+              HBondDonor* donor = new HBondDonor();
+              donor->donorAtom = atom1;
+              donor->donatedHydrogen = atom2;
+              hBondDonors_.push_back( donor );
+            }
+          }
+        }
+      }
+    }
+  }
+  
   RealType Molecule::getMass() {
     StuntDouble* sd;
     std::vector<StuntDouble*>::iterator i;
@@ -215,7 +302,26 @@ namespace OpenMD {
 
     return com;
   }
+  
+  Vector3d Molecule::getCom(int snapshotNo) {
+    StuntDouble* sd;
+    std::vector<StuntDouble*>::iterator i;
+    Vector3d com;
+    RealType totalMass = 0;
+    RealType mass;
+    
+    for (sd = beginIntegrableObject(i); sd != NULL; sd = 
+           nextIntegrableObject(i)){
+      mass = sd->getMass();
+      totalMass += mass;
+      com += sd->getPos(snapshotNo) * mass;    
+    }
+    
+    com /= totalMass;
 
+    return com;
+  }
+
   void Molecule::moveCom(const Vector3d& delta) {
     StuntDouble* sd;
     std::vector<StuntDouble*>::iterator i;