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 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
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, 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   */
42  
43   /**
# Line 52 | Line 53
53   #include "primitives/Molecule.hpp"
54   #include "utils/MemoryUtils.hpp"
55   #include "utils/simError.h"
56 + #include "utils/StringUtils.hpp"
57  
58 < namespace oopse {
59 <  Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
60 <    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
61 <
62 <    }
63 <
58 > namespace OpenMD {
59 >  Molecule::Molecule(int stampId, int globalIndex, const std::string& molName,
60 >                     int region) :
61 >                                   globalIndex_(globalIndex),
62 >                                   stampId_(stampId),
63 >                                   region_(region),
64 >                                   moleculeName_(molName),
65 >                                   constrainTotalCharge_(false) {
66 >  }
67 >  
68    Molecule::~Molecule() {
69 <
69 >    
70      MemoryUtils::deletePointers(atoms_);
71      MemoryUtils::deletePointers(bonds_);
72      MemoryUtils::deletePointers(bends_);
73      MemoryUtils::deletePointers(torsions_);
74 +    MemoryUtils::deletePointers(inversions_);
75      MemoryUtils::deletePointers(rigidBodies_);
76      MemoryUtils::deletePointers(cutoffGroups_);
77      MemoryUtils::deletePointers(constraintPairs_);
78      MemoryUtils::deletePointers(constraintElems_);
79 <    //integrableObjects_ don't own the objects
79 >
80 >    // integrableObjects_ don't own the objects
81      integrableObjects_.clear();
82 +    fluctuatingCharges_.clear();
83      
84    }
85 <
85 >  
86    void Molecule::addAtom(Atom* atom) {
87      if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
88        atoms_.push_back(atom);
89      }
90    }
91 <
91 >  
92    void Molecule::addBond(Bond* bond) {
93      if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
94        bonds_.push_back(bond);
95      }
96    }
97 <
97 >  
98    void Molecule::addBend(Bend* bend) {
99      if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
100        bends_.push_back(bend);
101      }
102    }
103 <
103 >  
104    void Molecule::addTorsion(Torsion* torsion) {
105 <    if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
105 >    if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
106 >        torsions_.end()) {
107        torsions_.push_back(torsion);
108      }
109    }
110  
111 +  void Molecule::addInversion(Inversion* inversion) {
112 +    if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
113 +        inversions_.end()) {
114 +      inversions_.push_back(inversion);
115 +    }
116 +  }
117 +  
118    void Molecule::addRigidBody(RigidBody *rb) {
119 <    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
119 >    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
120 >        rigidBodies_.end()) {
121        rigidBodies_.push_back(rb);
122      }
123    }
124 <
124 >  
125    void Molecule::addCutoffGroup(CutoffGroup* cp) {
126 <    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
126 >    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
127 >        cutoffGroups_.end()) {
128        cutoffGroups_.push_back(cp);
129 <    }
111 <
129 >    }    
130    }
131 <
131 >  
132    void Molecule::addConstraintPair(ConstraintPair* cp) {
133 <    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
133 >    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
134 >        constraintPairs_.end()) {
135        constraintPairs_.push_back(cp);
136 <    }
118 <
136 >    }    
137    }
138 <
138 >  
139    void Molecule::addConstraintElem(ConstraintElem* cp) {
140 <    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) {
140 >    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
141 >        constraintElems_.end()) {
142        constraintElems_.push_back(cp);
143      }
125
144    }
145 <
145 >  
146    void Molecule::complete() {
147      
148      std::set<Atom*> rigidAtoms;
149 +    Atom* atom;
150 +    Atom* atom1;
151 +    Atom* atom2;
152 +    AtomIterator ai, aj;
153      RigidBody* rb;
154 <    std::vector<RigidBody*>::iterator rbIter;
154 >    RigidBodyIterator rbIter;
155 >    Bond* bond;
156 >    BondIterator bi;
157  
158 <    
158 >    // Get list of all the atoms that are part of rigid bodies
159 >
160      for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
161        rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
162      }
163 +    
164 +    // add any atom that wasn't part of a rigid body to the list of integrableObjects
165  
139    Atom* atom;
140    AtomIterator ai;
166      for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
167 <  
168 <      if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
169 <        //if an atom does not belong to a rigid body, it is an integrable object
170 <        integrableObjects_.push_back(*ai);
167 >      
168 >      if (rigidAtoms.find(atom) == rigidAtoms.end()) {
169 >
170 >        // If an atom does not belong to a rigid body, it is an
171 >        // integrable object
172 >
173 >        integrableObjects_.push_back(atom);
174        }
175      }
176 +    
177 +    // then add the rigid bodies themselves to the integrableObjects
178  
149    //find all free atoms (which do not belong to rigid bodies)  
150    //performs the "difference" operation from set theory,  the output range contains a copy of every
151    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
152    //[rigidAtoms.begin(), rigidAtoms.end()).
153    //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
154    //                        std::back_inserter(integrableObjects_));
155
156    //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
157    //    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
158    //    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
159    //
160    //    painCave.isFatal = 1;
161    //    simError();        
162    //}
179      for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
180        integrableObjects_.push_back(rb);
181      }
166    //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
167  }
182  
183 +    // find the atoms that are fluctuating charges and add them to the
184 +    // fluctuatingCharges_ vector
185 +
186 +    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
187 +      if ( atom->isFluctuatingCharge() )
188 +        fluctuatingCharges_.push_back( atom );      
189 +    }
190 +
191 +
192 +    // find the electronegative atoms and add them to the
193 +    // hBondAcceptors_ vector:
194 +    
195 +    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
196 +      AtomType* at = atom->getAtomType();
197 +      // get the chain of base types for this atom type:
198 +      std::vector<AtomType*> ayb = at->allYourBase();
199 +      // use the last type in the chain of base types for the name:
200 +      std::string bn = UpperCase(ayb[ayb.size()-1]->getName());
201 +
202 +        if (bn.compare("O")==0 || bn.compare("N")==0
203 +            || bn.compare("F")==0)
204 +          hBondAcceptors_.push_back( atom );
205 +      
206 +    }
207 +    
208 +    // find electronegative atoms that are either bonded to
209 +    // hydrogens or are present in the same rigid bodies:
210 +    
211 +    for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) {
212 +      Atom* atom1 = bond->getAtomA();
213 +      Atom* atom2 = bond->getAtomB();
214 +      AtomType* at1 = atom1->getAtomType();
215 +      AtomType* at2 = atom1->getAtomType();
216 +      // get the chain of base types for this atom type:
217 +      std::vector<AtomType*> ayb1 = at1->allYourBase();
218 +      std::vector<AtomType*> ayb2 = at2->allYourBase();
219 +      // use the last type in the chain of base types for the name:
220 +      std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
221 +      std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
222 +      
223 +      if (bn1.compare("H")==0) {
224 +        if (bn2.compare("O")==0 || bn2.compare("N")==0
225 +            || bn2.compare("F")==0) {
226 +          HBondDonor* donor = new HBondDonor();
227 +          donor->donorAtom = atom2;
228 +          donor->donatedHydrogen = atom1;
229 +          hBondDonors_.push_back( donor );
230 +        }
231 +      }
232 +      if (bn2.compare("H")==0) {
233 +        if (bn1.compare("O")==0 || bn1.compare("N")==0
234 +            || bn1.compare("F")==0) {
235 +          HBondDonor* donor = new HBondDonor();
236 +          donor->donorAtom = atom1;
237 +          donor->donatedHydrogen = atom2;
238 +            hBondDonors_.push_back( donor );
239 +        }
240 +      }
241 +    }
242 +    
243 +    for (rb = beginRigidBody(rbIter); rb != NULL;
244 +         rb = nextRigidBody(rbIter)) {
245 +      for(atom1 = rb->beginAtom(ai); atom1 != NULL;
246 +          atom1 = rb->nextAtom(ai)) {
247 +        AtomType* at1 = atom1->getAtomType();
248 +        // get the chain of base types for this atom type:
249 +        std::vector<AtomType*> ayb1 = at1->allYourBase();
250 +        // use the last type in the chain of base types for the name:
251 +        std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
252 +        
253 +        if (bn1.compare("O")==0 || bn1.compare("N")==0
254 +            || bn1.compare("F")==0) {
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 = UpperCase(ayb2[ayb2.size()-1]->getName());
262 +            if (bn2.compare("H")==0) {              
263 +              HBondDonor* donor = new HBondDonor();
264 +              donor->donorAtom = atom1;
265 +              donor->donatedHydrogen = atom2;
266 +              hBondDonors_.push_back( donor );
267 +            }
268 +          }
269 +        }
270 +      }
271 +    }
272 +  }
273 +  
274    RealType Molecule::getMass() {
275      StuntDouble* sd;
276      std::vector<StuntDouble*>::iterator i;
277      RealType mass = 0.0;
278 <
279 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
278 >    
279 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
280 >           nextIntegrableObject(i)){
281        mass += sd->getMass();
282      }
283 <
284 <    return mass;
179 <
283 >    
284 >    return mass;    
285    }
286  
287    Vector3d Molecule::getCom() {
# Line 186 | Line 291 | namespace oopse {
291      RealType totalMass = 0;
292      RealType mass;
293      
294 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
294 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
295 >           nextIntegrableObject(i)){
296        mass = sd->getMass();
297        totalMass += mass;
298        com += sd->getPos() * mass;    
299      }
300 +    
301 +    com /= totalMass;
302  
303 +    return com;
304 +  }
305 +  
306 +  Vector3d Molecule::getCom(int snapshotNo) {
307 +    StuntDouble* sd;
308 +    std::vector<StuntDouble*>::iterator i;
309 +    Vector3d com;
310 +    RealType totalMass = 0;
311 +    RealType mass;
312 +    
313 +    for (sd = beginIntegrableObject(i); sd != NULL; sd =
314 +           nextIntegrableObject(i)){
315 +      mass = sd->getMass();
316 +      totalMass += mass;
317 +      com += sd->getPos(snapshotNo) * mass;    
318 +    }
319 +    
320      com /= totalMass;
321  
322      return com;
# Line 201 | Line 326 | namespace oopse {
326      StuntDouble* sd;
327      std::vector<StuntDouble*>::iterator i;
328      
329 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
329 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
330 >           nextIntegrableObject(i)){
331        sd->setPos(sd->getPos() + delta);
332 <    }
207 <
332 >    }    
333    }
334  
335    Vector3d Molecule::getComVel() {
# Line 214 | Line 339 | namespace oopse {
339      RealType totalMass = 0;
340      RealType mass;
341      
342 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
342 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
343 >           nextIntegrableObject(i)){
344        mass = sd->getMass();
345        totalMass += mass;
346        velCom += sd->getVel() * mass;    
347      }
348 <
348 >    
349      velCom /= totalMass;
350 <
350 >    
351      return velCom;
352    }
353  
# Line 230 | Line 356 | namespace oopse {
356      Bond* bond;
357      Bend* bend;
358      Torsion* torsion;
359 +    Inversion* inversion;
360      Molecule::BondIterator bondIter;;
361      Molecule::BendIterator  bendIter;
362      Molecule::TorsionIterator  torsionIter;
363 +    Molecule::InversionIterator  inversionIter;
364  
365      RealType potential = 0.0;
366  
# Line 244 | Line 372 | namespace oopse {
372        potential += bend->getPotential();
373      }
374  
375 <    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
375 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
376 >           nextTorsion(torsionIter)) {
377        potential += torsion->getPotential();
378      }
379 <
379 >    
380 >    for (inversion = beginInversion(inversionIter); torsion != NULL;
381 >         inversion =  nextInversion(inversionIter)) {
382 >      potential += inversion->getPotential();
383 >    }
384 >    
385      return potential;
386 +    
387 +  }
388 +  
389 +  void Molecule::addProperty(GenericData* genData) {
390 +    properties_.addProperty(genData);  
391 +  }
392  
393 +  void Molecule::removeProperty(const std::string& propName) {
394 +    properties_.removeProperty(propName);  
395    }
396  
397 +  void Molecule::clearProperties() {
398 +    properties_.clearProperties();
399 +  }
400 +
401 +  std::vector<std::string> Molecule::getPropertyNames() {
402 +    return properties_.getPropertyNames();  
403 +  }
404 +      
405 +  std::vector<GenericData*> Molecule::getProperties() {
406 +    return properties_.getProperties();
407 +  }
408 +
409 +  GenericData* Molecule::getPropertyByName(const std::string& propName) {
410 +    return properties_.getPropertyByName(propName);
411 +  }
412 +
413    std::ostream& operator <<(std::ostream& o, Molecule& mol) {
414      o << std::endl;
415      o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
# Line 259 | Line 417 | namespace oopse {
417      o << mol.getNBonds() << " bonds" << std::endl;
418      o << mol.getNBends() << " bends" << std::endl;
419      o << mol.getNTorsions() << " torsions" << std::endl;
420 +    o << mol.getNInversions() << " inversions" << std::endl;
421      o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
422 <    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
423 <    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
424 <    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
422 >    o << mol.getNIntegrableObjects() << " integrable objects" << std::endl;
423 >    o << mol.getNCutoffGroups() << " cutoff groups" << std::endl;
424 >    o << mol.getNConstraintPairs() << " constraint pairs" << std::endl;
425 >    o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl;
426      return o;
427    }
428 <
429 < }//end namespace oopse
428 >  
429 > }//end namespace OpenMD

Comparing trunk/src/primitives/Molecule.cpp (property svn:keywords):
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines