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 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2053 by gezelter, Mon Jan 12 16:04:43 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/ElementsTable.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) : stampId_(stampId),
61 >                                   globalIndex_(globalIndex),
62 >                                   moleculeName_(molName),
63 >                                   region_(region),
64 >                                   constrainTotalCharge_(false) {
65 >  }
66 >  
67    Molecule::~Molecule() {
68 <
68 >    
69      MemoryUtils::deletePointers(atoms_);
70      MemoryUtils::deletePointers(bonds_);
71      MemoryUtils::deletePointers(bends_);
72      MemoryUtils::deletePointers(torsions_);
73 +    MemoryUtils::deletePointers(inversions_);
74      MemoryUtils::deletePointers(rigidBodies_);
75      MemoryUtils::deletePointers(cutoffGroups_);
76      MemoryUtils::deletePointers(constraintPairs_);
77      MemoryUtils::deletePointers(constraintElems_);
78 <    //integrableObjects_ don't own the objects
78 >
79 >    // integrableObjects_ don't own the objects
80      integrableObjects_.clear();
81 +    fluctuatingCharges_.clear();
82      
83    }
84 <
84 >  
85    void Molecule::addAtom(Atom* atom) {
86      if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
87        atoms_.push_back(atom);
88      }
89    }
90 <
90 >  
91    void Molecule::addBond(Bond* bond) {
92      if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
93        bonds_.push_back(bond);
94      }
95    }
96 <
96 >  
97    void Molecule::addBend(Bend* bend) {
98      if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
99        bends_.push_back(bend);
100      }
101    }
102 <
102 >  
103    void Molecule::addTorsion(Torsion* torsion) {
104 <    if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
104 >    if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
105 >        torsions_.end()) {
106        torsions_.push_back(torsion);
107      }
108    }
109  
110 +  void Molecule::addInversion(Inversion* inversion) {
111 +    if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
112 +        inversions_.end()) {
113 +      inversions_.push_back(inversion);
114 +    }
115 +  }
116 +  
117    void Molecule::addRigidBody(RigidBody *rb) {
118 <    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
118 >    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
119 >        rigidBodies_.end()) {
120        rigidBodies_.push_back(rb);
121      }
122    }
123 <
123 >  
124    void Molecule::addCutoffGroup(CutoffGroup* cp) {
125 <    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
125 >    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
126 >        cutoffGroups_.end()) {
127        cutoffGroups_.push_back(cp);
128 <    }
111 <
128 >    }    
129    }
130 <
130 >  
131    void Molecule::addConstraintPair(ConstraintPair* cp) {
132 <    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
132 >    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
133 >        constraintPairs_.end()) {
134        constraintPairs_.push_back(cp);
135 <    }
118 <
135 >    }    
136    }
137 <
137 >  
138    void Molecule::addConstraintElem(ConstraintElem* cp) {
139 <    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) {
139 >    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
140 >        constraintElems_.end()) {
141        constraintElems_.push_back(cp);
142      }
125
143    }
144 <
144 >  
145    void Molecule::complete() {
146      
147      std::set<Atom*> rigidAtoms;
148 +    Atom* atom;
149 +    Atom* atom1;
150 +    Atom* atom2;
151 +    AtomIterator ai, aj;
152      RigidBody* rb;
153 <    std::vector<RigidBody*>::iterator rbIter;
153 >    RigidBodyIterator rbIter;
154 >    Bond* bond;
155 >    BondIterator bi;
156  
157 <    
157 >    // Get list of all the atoms that are part of rigid bodies
158 >
159      for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
160        rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
161      }
162 +    
163 +    // add any atom that wasn't part of a rigid body to the list of integrableObjects
164  
139    Atom* atom;
140    AtomIterator ai;
165      for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
166 <  
167 <      if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
168 <        //if an atom does not belong to a rigid body, it is an integrable object
169 <        integrableObjects_.push_back(*ai);
166 >      
167 >      if (rigidAtoms.find(atom) == rigidAtoms.end()) {
168 >
169 >        // If an atom does not belong to a rigid body, it is an
170 >        // integrable object
171 >
172 >        integrableObjects_.push_back(atom);
173        }
174      }
175 +    
176 +    // then add the rigid bodies themselves to the integrableObjects
177  
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    //}
178      for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
179        integrableObjects_.push_back(rb);
180      }
166    //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
167  }
181  
182 <  double Molecule::getMass() {
183 <    StuntDouble* sd;
171 <    std::vector<StuntDouble*>::iterator i;
172 <    double mass = 0.0;
182 >    // find the atoms that are fluctuating charges and add them to the
183 >    // fluctuatingCharges_ vector
184  
185 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
186 <      mass += sd->getMass();
185 >    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
186 >      if ( atom->isFluctuatingCharge() )
187 >        fluctuatingCharges_.push_back( atom );      
188      }
189  
178    return mass;
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;
275 +    RealType mass = 0.0;
276 +    
277 +    for (sd = beginIntegrableObject(i); sd != NULL; sd =
278 +           nextIntegrableObject(i)){
279 +      mass += sd->getMass();
280 +    }
281 +    
282 +    return mass;    
283 +  }
284  
285    Vector3d Molecule::getCom() {
286      StuntDouble* sd;
287      std::vector<StuntDouble*>::iterator i;
288      Vector3d com;
289 <    double totalMass = 0;
290 <    double mass;
289 >    RealType totalMass = 0;
290 >    RealType mass;
291      
292 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
292 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
293 >           nextIntegrableObject(i)){
294        mass = sd->getMass();
295        totalMass += mass;
296        com += sd->getPos() * mass;    
297      }
298 +    
299 +    com /= totalMass;
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;
# Line 201 | Line 324 | namespace oopse {
324      StuntDouble* sd;
325      std::vector<StuntDouble*>::iterator i;
326      
327 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
327 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
328 >           nextIntegrableObject(i)){
329        sd->setPos(sd->getPos() + delta);
330 <    }
207 <
330 >    }    
331    }
332  
333    Vector3d Molecule::getComVel() {
334      StuntDouble* sd;
335      std::vector<StuntDouble*>::iterator i;
336      Vector3d velCom;
337 <    double totalMass = 0;
338 <    double mass;
337 >    RealType totalMass = 0;
338 >    RealType mass;
339      
340 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
340 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
341 >           nextIntegrableObject(i)){
342        mass = sd->getMass();
343        totalMass += mass;
344        velCom += sd->getVel() * mass;    
345      }
346 <
346 >    
347      velCom /= totalMass;
348 <
348 >    
349      return velCom;
350    }
351  
352 <  double Molecule::getPotential() {
352 >  RealType Molecule::getPotential() {
353  
354      Bond* bond;
355      Bend* bend;
356      Torsion* torsion;
357 +    Inversion* inversion;
358      Molecule::BondIterator bondIter;;
359      Molecule::BendIterator  bendIter;
360      Molecule::TorsionIterator  torsionIter;
361 +    Molecule::InversionIterator  inversionIter;
362  
363 <    double potential = 0.0;
363 >    RealType potential = 0.0;
364  
365      for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
366        potential += bond->getPotential();
# Line 244 | Line 370 | namespace oopse {
370        potential += bend->getPotential();
371      }
372  
373 <    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
373 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
374 >           nextTorsion(torsionIter)) {
375        potential += torsion->getPotential();
376      }
377 <
377 >    
378 >    for (inversion = beginInversion(inversionIter); torsion != NULL;
379 >         inversion =  nextInversion(inversionIter)) {
380 >      potential += inversion->getPotential();
381 >    }
382 >    
383      return potential;
384 +    
385 +  }
386 +  
387 +  void Molecule::addProperty(GenericData* genData) {
388 +    properties_.addProperty(genData);  
389 +  }
390  
391 +  void Molecule::removeProperty(const std::string& propName) {
392 +    properties_.removeProperty(propName);  
393    }
394  
395 +  void Molecule::clearProperties() {
396 +    properties_.clearProperties();
397 +  }
398 +
399 +  std::vector<std::string> Molecule::getPropertyNames() {
400 +    return properties_.getPropertyNames();  
401 +  }
402 +      
403 +  std::vector<GenericData*> Molecule::getProperties() {
404 +    return properties_.getProperties();
405 +  }
406 +
407 +  GenericData* Molecule::getPropertyByName(const std::string& propName) {
408 +    return properties_.getPropertyByName(propName);
409 +  }
410 +
411    std::ostream& operator <<(std::ostream& o, Molecule& mol) {
412      o << std::endl;
413      o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
# Line 259 | Line 415 | namespace oopse {
415      o << mol.getNBonds() << " bonds" << std::endl;
416      o << mol.getNBends() << " bends" << std::endl;
417      o << mol.getNTorsions() << " torsions" << std::endl;
418 +    o << mol.getNInversions() << " inversions" << std::endl;
419      o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
420 <    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
421 <    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
422 <    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
420 >    o << mol.getNIntegrableObjects() << " integrable objects" << std::endl;
421 >    o << mol.getNCutoffGroups() << " cutoff groups" << std::endl;
422 >    o << mol.getNConstraintPairs() << " constraint pairs" << std::endl;
423 >    o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl;
424      return o;
425    }
426 <
427 < }//end namespace oopse
426 >  
427 > }//end namespace OpenMD

Comparing trunk/src/primitives/Molecule.cpp (property svn:keywords):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2053 by gezelter, Mon Jan 12 16:04:43 2015 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines