ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/primitives/Molecule.cpp
(Generate patch)

Comparing trunk/src/primitives/Molecule.cpp (file contents):
Revision 2052 by gezelter, Fri Jan 9 19:06:35 2015 UTC vs.
Revision 2054 by gezelter, Tue Jan 13 20:14:57 2015 UTC

# Line 53 | Line 53
53   #include "primitives/Molecule.hpp"
54   #include "utils/MemoryUtils.hpp"
55   #include "utils/simError.h"
56 < #include "utils/ElementsTable.hpp"
56 > #include "utils/StringUtils.hpp"
57  
58   namespace OpenMD {
59    Molecule::Molecule(int stampId, int globalIndex, const std::string& molName,
# Line 196 | Line 196 | namespace OpenMD {
196        // get the chain of base types for this atom type:
197        std::vector<AtomType*> ayb = at->allYourBase();
198        // use the last type in the chain of base types for the name:
199 <      std::string bn = ayb[ayb.size()-1]->getName();
199 >      std::string bn = UpperCase(ayb[ayb.size()-1]->getName());
200  
201 <      int obanum = etab.GetAtomicNum(bn.c_str());
202 <      if (obanum != 0) {
203 <        RealType eneg = etab.GetElectroNeg(obanum);
204 <        if (eneg > 3.01) {
201 >        if (bn.compare("O")==0 || bn.compare("N")==0
202 >            || bn.compare("F")==0)
203            hBondAcceptors_.push_back( atom );
204 <        }
207 <      }
204 >      
205      }
209
210    // find electronegative atoms that are either bonded to hydrogens or are
211    // present in the same rigid bodies:
206      
207 +    // find electronegative atoms that are either bonded to
208 +    // hydrogens or are present in the same rigid bodies:
209 +    
210      for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) {
211        Atom* atom1 = bond->getAtomA();
212        Atom* atom2 = bond->getAtomB();
# Line 219 | Line 216 | namespace OpenMD {
216        std::vector<AtomType*> ayb1 = at1->allYourBase();
217        std::vector<AtomType*> ayb2 = at2->allYourBase();
218        // use the last type in the chain of base types for the name:
219 <      std::string bn1 = ayb1[ayb1.size()-1]->getName();
220 <      std::string bn2 = ayb2[ayb2.size()-1]->getName();
221 <      int obanum1 = etab.GetAtomicNum(bn1.c_str());
222 <      int obanum2 = etab.GetAtomicNum(bn2.c_str());
223 <
224 <      if (obanum1 == 1) {              
225 <        if (obanum2 != 0) {
226 <          RealType eneg = etab.GetElectroNeg(obanum2);
227 <          if (eneg > 3.01) {
228 <            HBondDonor* donor = new HBondDonor();
232 <            donor->donorAtom = atom2;
233 <            donor->donatedHydrogen = atom1;
234 <            hBondDonors_.push_back( donor );
235 <          }
219 >      std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
220 >      std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
221 >      
222 >      if (bn1.compare("H")==0) {
223 >        if (bn2.compare("O")==0 || bn2.compare("N")==0
224 >            || bn2.compare("F")==0) {
225 >          HBondDonor* donor = new HBondDonor();
226 >          donor->donorAtom = atom2;
227 >          donor->donatedHydrogen = atom1;
228 >          hBondDonors_.push_back( donor );
229          }
230        }
231 <      if (obanum2 == 1) {
232 <        if (obanum1 != 0) {
233 <          RealType eneg = etab.GetElectroNeg(obanum1);
234 <          if (eneg > 3.01) {
235 <            HBondDonor* donor = new HBondDonor();
236 <            donor->donorAtom = atom1;
237 <            donor->donatedHydrogen = atom2;
238 <            hBondDonors_.push_back( donor );
239 <          }
247 <        }
248 <      }
231 >      if (bn2.compare("H")==0) {
232 >        if (bn1.compare("O")==0 || bn1.compare("N")==0
233 >            || bn1.compare("F")==0) {
234 >          HBondDonor* donor = new HBondDonor();
235 >          donor->donorAtom = atom1;
236 >          donor->donatedHydrogen = atom2;
237 >            hBondDonors_.push_back( donor );
238 >        }
239 >      }
240      }
241 <
242 <    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
243 <      for(atom1 = rb->beginAtom(ai); atom1 != NULL; atom1 = rb->nextAtom(ai)) {
241 >    
242 >    for (rb = beginRigidBody(rbIter); rb != NULL;
243 >         rb = nextRigidBody(rbIter)) {
244 >      for(atom1 = rb->beginAtom(ai); atom1 != NULL;
245 >          atom1 = rb->nextAtom(ai)) {
246          AtomType* at1 = atom1->getAtomType();
247          // get the chain of base types for this atom type:
248          std::vector<AtomType*> ayb1 = at1->allYourBase();
249          // use the last type in the chain of base types for the name:
250 <        std::string bn1 = ayb1[ayb1.size()-1]->getName();
251 <        int obanum1 = etab.GetAtomicNum(bn1.c_str());
252 <        if (obanum1 != 0) {
253 <          RealType eneg = etab.GetElectroNeg(obanum1);
254 <          if (eneg > 3.01) {
255 <            for(atom2 = rb->beginAtom(aj); atom2 != NULL;
256 <                atom2 = rb->nextAtom(aj)) {
257 <              AtomType* at2 = atom2->getAtomType();
258 <              // get the chain of base types for this atom type:              
259 <              std::vector<AtomType*> ayb2 = at2->allYourBase();
260 <              // use the last type in the chain of base types for the name:
261 <              std::string bn2 = ayb2[ayb2.size()-1]->getName();
262 <              int obanum2 = etab.GetAtomicNum(bn2.c_str());
263 <              if (obanum2 == 1) {
264 <                HBondDonor* donor = new HBondDonor();
265 <                donor->donorAtom = atom1;
273 <                donor->donatedHydrogen = atom2;
274 <                hBondDonors_.push_back( donor );
275 <              }
250 >        std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
251 >        
252 >        if (bn1.compare("O")==0 || bn1.compare("N")==0
253 >            || bn1.compare("F")==0) {
254 >          for(atom2 = rb->beginAtom(aj); atom2 != NULL;
255 >              atom2 = rb->nextAtom(aj)) {
256 >            AtomType* at2 = atom2->getAtomType();
257 >            // get the chain of base types for this atom type:              
258 >            std::vector<AtomType*> ayb2 = at2->allYourBase();
259 >            // use the last type in the chain of base types for the name:
260 >            std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
261 >            if (bn2.compare("H")==0) {              
262 >              HBondDonor* donor = new HBondDonor();
263 >              donor->donorAtom = atom1;
264 >              donor->donatedHydrogen = atom2;
265 >              hBondDonors_.push_back( donor );
266              }
267            }
268          }
269        }
270 <    }          
270 >    }
271    }
272 <
272 >  
273    RealType Molecule::getMass() {
274      StuntDouble* sd;
275      std::vector<StuntDouble*>::iterator i;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines