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 1782 by gezelter, Wed Aug 22 02:28:28 2012 UTC vs.
Revision 2053 by gezelter, Mon Jan 12 16:04:43 2015 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39   * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40   * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
# 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 179 | Line 187 | namespace OpenMD {
187          fluctuatingCharges_.push_back( atom );      
188      }
189  
182 –  }
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 +      if (bn.compare("O")==0 || bn.compare("N")==0 || bn.compare("F")==0)
202 +        hBondAcceptors_.push_back( atom );
203 +
204 +    }
205 +    
206 +    // find electronegative atoms that are either bonded to
207 +    // hydrogens or are present in the same rigid bodies:
208 +    
209 +    for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) {
210 +      Atom* atom1 = bond->getAtomA();
211 +      Atom* atom2 = bond->getAtomB();
212 +      AtomType* at1 = atom1->getAtomType();
213 +      AtomType* at2 = atom1->getAtomType();
214 +      // get the chain of base types for this atom type:
215 +      std::vector<AtomType*> ayb1 = at1->allYourBase();
216 +      std::vector<AtomType*> ayb2 = at2->allYourBase();
217 +      // use the last type in the chain of base types for the name:
218 +      std::string bn1 = ayb1[ayb1.size()-1]->getName();
219 +      std::string bn2 = ayb2[ayb2.size()-1]->getName();
220 +      
221 +      if (bn1.compare("H")==0) {
222 +        if (bn2.compare("O")==0 || bn2.compare("N")==0
223 +            || bn2.compare("F")==0) {
224 +          HBondDonor* donor = new HBondDonor();
225 +          donor->donorAtom = atom2;
226 +          donor->donatedHydrogen = atom1;
227 +          hBondDonors_.push_back( donor );
228 +        }
229 +      }
230 +      if (bn2.compare("H")==0) {
231 +        if (bn1.compare("O")==0 || bn1.compare("N")==0
232 +            || bn1.compare("F")==0) {
233 +          HBondDonor* donor = new HBondDonor();
234 +          donor->donorAtom = atom1;
235 +          donor->donatedHydrogen = atom2;
236 +            hBondDonors_.push_back( donor );
237 +        }
238 +      }
239 +    }
240 +    
241 +    for (rb = beginRigidBody(rbIter); rb != NULL;
242 +         rb = nextRigidBody(rbIter)) {
243 +      for(atom1 = rb->beginAtom(ai); atom1 != NULL;
244 +          atom1 = rb->nextAtom(ai)) {
245 +        AtomType* at1 = atom1->getAtomType();
246 +        // get the chain of base types for this atom type:
247 +        std::vector<AtomType*> ayb1 = at1->allYourBase();
248 +        // use the last type in the chain of base types for the name:
249 +        std::string bn1 = ayb1[ayb1.size()-1]->getName();
250 +        
251 +        if (bn1.compare("O")==0 || bn1.compare("N")==0
252 +            || bn1.compare("F")==0) {
253 +          for(atom2 = rb->beginAtom(aj); atom2 != NULL;
254 +              atom2 = rb->nextAtom(aj)) {
255 +            AtomType* at2 = atom2->getAtomType();
256 +            // get the chain of base types for this atom type:              
257 +            std::vector<AtomType*> ayb2 = at2->allYourBase();
258 +            // use the last type in the chain of base types for the name:
259 +            std::string bn2 = ayb2[ayb2.size()-1]->getName();
260 +            if (bn2.compare("H")==0) {              
261 +              HBondDonor* donor = new HBondDonor();
262 +              donor->donorAtom = atom1;
263 +              donor->donatedHydrogen = atom2;
264 +              hBondDonors_.push_back( donor );
265 +            }
266 +          }
267 +        }
268 +      }
269 +    }
270 +  }
271 +  
272    RealType Molecule::getMass() {
273      StuntDouble* sd;
274      std::vector<StuntDouble*>::iterator i;
# Line 212 | Line 300 | namespace OpenMD {
300  
301      return com;
302    }
303 +  
304 +  Vector3d Molecule::getCom(int snapshotNo) {
305 +    StuntDouble* sd;
306 +    std::vector<StuntDouble*>::iterator i;
307 +    Vector3d com;
308 +    RealType totalMass = 0;
309 +    RealType mass;
310 +    
311 +    for (sd = beginIntegrableObject(i); sd != NULL; sd =
312 +           nextIntegrableObject(i)){
313 +      mass = sd->getMass();
314 +      totalMass += mass;
315 +      com += sd->getPos(snapshotNo) * mass;    
316 +    }
317 +    
318 +    com /= totalMass;
319  
320 +    return com;
321 +  }
322 +
323    void Molecule::moveCom(const Vector3d& delta) {
324      StuntDouble* sd;
325      std::vector<StuntDouble*>::iterator i;

Diff Legend

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