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 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2052 by gezelter, Fri Jan 9 19:06:35 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"
57  
58   namespace OpenMD {
59 <  Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
60 <    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName),
61 <      constrainTotalCharge_(false) {
59 >  Molecule::Molecule(int stampId, int globalIndex, const std::string& molName,
60 >                     int region) : stampId_(stampId),
61 >                                   globalIndex_(globalIndex),
62 >                                   moleculeName_(molName),
63 >                                   region_(region),
64 >                                   constrainTotalCharge_(false) {
65    }
66    
67    Molecule::~Molecule() {
# Line 142 | Line 146 | namespace OpenMD {
146      
147      std::set<Atom*> rigidAtoms;
148      Atom* atom;
149 <    AtomIterator ai;
149 >    Atom* atom1;
150 >    Atom* atom2;
151 >    AtomIterator ai, aj;
152      RigidBody* rb;
153      RigidBodyIterator rbIter;
154 +    Bond* bond;
155 +    BondIterator bi;
156  
157      // Get list of all the atoms that are part of rigid bodies
158  
# Line 177 | Line 185 | namespace OpenMD {
185      for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
186        if ( atom->isFluctuatingCharge() )
187          fluctuatingCharges_.push_back( atom );      
188 +    }
189 +
190 +
191 +    // find the electronegative atoms and add them to the
192 +    // hBondAcceptors_ vector:
193 +    
194 +    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
195 +      AtomType* at = atom->getAtomType();
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();
200 +
201 +      int obanum = etab.GetAtomicNum(bn.c_str());
202 +      if (obanum != 0) {
203 +        RealType eneg = etab.GetElectroNeg(obanum);
204 +        if (eneg > 3.01) {
205 +          hBondAcceptors_.push_back( atom );
206 +        }
207 +      }
208 +    }
209 +
210 +    // find electronegative atoms that are either bonded to hydrogens or are
211 +    // present in the same rigid bodies:
212 +    
213 +    for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) {
214 +      Atom* atom1 = bond->getAtomA();
215 +      Atom* atom2 = bond->getAtomB();
216 +      AtomType* at1 = atom1->getAtomType();
217 +      AtomType* at2 = atom1->getAtomType();
218 +      // get the chain of base types for this atom type:
219 +      std::vector<AtomType*> ayb1 = at1->allYourBase();
220 +      std::vector<AtomType*> ayb2 = at2->allYourBase();
221 +      // use the last type in the chain of base types for the name:
222 +      std::string bn1 = ayb1[ayb1.size()-1]->getName();
223 +      std::string bn2 = ayb2[ayb2.size()-1]->getName();
224 +      int obanum1 = etab.GetAtomicNum(bn1.c_str());
225 +      int obanum2 = etab.GetAtomicNum(bn2.c_str());
226 +
227 +      if (obanum1 == 1) {              
228 +        if (obanum2 != 0) {
229 +          RealType eneg = etab.GetElectroNeg(obanum2);
230 +          if (eneg > 3.01) {
231 +            HBondDonor* donor = new HBondDonor();
232 +            donor->donorAtom = atom2;
233 +            donor->donatedHydrogen = atom1;
234 +            hBondDonors_.push_back( donor );
235 +          }
236 +        }
237 +      }
238 +      if (obanum2 == 1) {
239 +        if (obanum1 != 0) {
240 +          RealType eneg = etab.GetElectroNeg(obanum1);
241 +          if (eneg > 3.01) {
242 +            HBondDonor* donor = new HBondDonor();
243 +            donor->donorAtom = atom1;
244 +            donor->donatedHydrogen = atom2;
245 +            hBondDonors_.push_back( donor );
246 +          }
247 +        }
248 +      }
249      }
250  
251 +    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
252 +      for(atom1 = rb->beginAtom(ai); atom1 != NULL; atom1 = rb->nextAtom(ai)) {
253 +        AtomType* at1 = atom1->getAtomType();
254 +        // get the chain of base types for this atom type:
255 +        std::vector<AtomType*> ayb1 = at1->allYourBase();
256 +        // use the last type in the chain of base types for the name:
257 +        std::string bn1 = ayb1[ayb1.size()-1]->getName();
258 +        int obanum1 = etab.GetAtomicNum(bn1.c_str());
259 +        if (obanum1 != 0) {
260 +          RealType eneg = etab.GetElectroNeg(obanum1);
261 +          if (eneg > 3.01) {
262 +            for(atom2 = rb->beginAtom(aj); atom2 != NULL;
263 +                atom2 = rb->nextAtom(aj)) {
264 +              AtomType* at2 = atom2->getAtomType();
265 +              // get the chain of base types for this atom type:              
266 +              std::vector<AtomType*> ayb2 = at2->allYourBase();
267 +              // use the last type in the chain of base types for the name:
268 +              std::string bn2 = ayb2[ayb2.size()-1]->getName();
269 +              int obanum2 = etab.GetAtomicNum(bn2.c_str());
270 +              if (obanum2 == 1) {
271 +                HBondDonor* donor = new HBondDonor();
272 +                donor->donorAtom = atom1;
273 +                donor->donatedHydrogen = atom2;
274 +                hBondDonors_.push_back( donor );
275 +              }
276 +            }
277 +          }
278 +        }
279 +      }
280 +    }          
281    }
282  
283    RealType Molecule::getMass() {
# Line 212 | Line 311 | namespace OpenMD {
311  
312      return com;
313    }
314 +  
315 +  Vector3d Molecule::getCom(int snapshotNo) {
316 +    StuntDouble* sd;
317 +    std::vector<StuntDouble*>::iterator i;
318 +    Vector3d com;
319 +    RealType totalMass = 0;
320 +    RealType mass;
321 +    
322 +    for (sd = beginIntegrableObject(i); sd != NULL; sd =
323 +           nextIntegrableObject(i)){
324 +      mass = sd->getMass();
325 +      totalMass += mass;
326 +      com += sd->getPos(snapshotNo) * mass;    
327 +    }
328 +    
329 +    com /= totalMass;
330  
331 +    return com;
332 +  }
333 +
334    void Molecule::moveCom(const Vector3d& delta) {
335      StuntDouble* sd;
336      std::vector<StuntDouble*>::iterator i;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines