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

Comparing:
trunk/src/primitives/Molecule.cpp (file contents), Revision 372 by tim, Mon Feb 21 15:28:56 2005 UTC vs.
branches/development/src/primitives/Molecule.cpp (file contents), Revision 1718 by gezelter, Thu May 24 01:29:59 2012 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# 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, 24107 (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 53 | Line 54
54   #include "utils/MemoryUtils.hpp"
55   #include "utils/simError.h"
56  
57 < namespace oopse {
58 < Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
59 <    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
57 > namespace OpenMD {
58 >  Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
59 >    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName),
60 >      constrainTotalCharge_(false) {
61 >  }
62 >  
63 >  Molecule::~Molecule() {
64 >    
65 >    MemoryUtils::deletePointers(atoms_);
66 >    MemoryUtils::deletePointers(bonds_);
67 >    MemoryUtils::deletePointers(bends_);
68 >    MemoryUtils::deletePointers(torsions_);
69 >    MemoryUtils::deletePointers(inversions_);
70 >    MemoryUtils::deletePointers(rigidBodies_);
71 >    MemoryUtils::deletePointers(cutoffGroups_);
72 >    MemoryUtils::deletePointers(constraintPairs_);
73 >    MemoryUtils::deletePointers(constraintElems_);
74  
75 < }
61 <
62 < Molecule::~Molecule() {
63 <
64 <    MemoryUtils::deleteVectorOfPointer(atoms_);
65 <    MemoryUtils::deleteVectorOfPointer(bonds_);
66 <    MemoryUtils::deleteVectorOfPointer(bends_);
67 <    MemoryUtils::deleteVectorOfPointer(torsions_);
68 <    MemoryUtils::deleteVectorOfPointer(rigidBodies_);
69 <    MemoryUtils::deleteVectorOfPointer(cutoffGroups_);
70 <    MemoryUtils::deleteVectorOfPointer(constraintPairs_);
71 <    MemoryUtils::deleteVectorOfPointer(constraintElems_);
72 <    //integrableObjects_ don't own the objects
75 >    // integrableObjects_ don't own the objects
76      integrableObjects_.clear();
77 +    fluctuatingCharges_.clear();
78      
79 < }
80 <
81 < void Molecule::addAtom(Atom* atom) {
79 >  }
80 >  
81 >  void Molecule::addAtom(Atom* atom) {
82      if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
83 <        atoms_.push_back(atom);
83 >      atoms_.push_back(atom);
84      }
85 < }
86 <
87 < void Molecule::addBond(Bond* bond) {
85 >  }
86 >  
87 >  void Molecule::addBond(Bond* bond) {
88      if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
89 <        bonds_.push_back(bond);
89 >      bonds_.push_back(bond);
90      }
91 < }
92 <
93 < void Molecule::addBend(Bend* bend) {
91 >  }
92 >  
93 >  void Molecule::addBend(Bend* bend) {
94      if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
95 <        bends_.push_back(bend);
95 >      bends_.push_back(bend);
96      }
97 < }
98 <
99 < void Molecule::addTorsion(Torsion* torsion) {
100 <    if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
101 <        torsions_.push_back(torsion);
97 >  }
98 >  
99 >  void Molecule::addTorsion(Torsion* torsion) {
100 >    if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
101 >        torsions_.end()) {
102 >      torsions_.push_back(torsion);
103      }
104 < }
104 >  }
105  
106 < void Molecule::addRigidBody(RigidBody *rb) {
107 <    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
108 <        rigidBodies_.push_back(rb);
106 >  void Molecule::addInversion(Inversion* inversion) {
107 >    if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
108 >        inversions_.end()) {
109 >      inversions_.push_back(inversion);
110      }
111 < }
112 <
113 < void Molecule::addCutoffGroup(CutoffGroup* cp) {
114 <    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
115 <        cutoffGroups_.push_back(cp);
111 >  }
112 >  
113 >  void Molecule::addRigidBody(RigidBody *rb) {
114 >    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
115 >        rigidBodies_.end()) {
116 >      rigidBodies_.push_back(rb);
117      }
118 <
119 < }
120 <
121 < void Molecule::addConstraintPair(ConstraintPair* cp) {
122 <    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
123 <        constraintPairs_.push_back(cp);
118 >  }
119 >  
120 >  void Molecule::addCutoffGroup(CutoffGroup* cp) {
121 >    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
122 >        cutoffGroups_.end()) {
123 >      cutoffGroups_.push_back(cp);
124 >    }    
125 >  }
126 >  
127 >  void Molecule::addConstraintPair(ConstraintPair* cp) {
128 >    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
129 >        constraintPairs_.end()) {
130 >      constraintPairs_.push_back(cp);
131 >    }    
132 >  }
133 >  
134 >  void Molecule::addConstraintElem(ConstraintElem* cp) {
135 >    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
136 >        constraintElems_.end()) {
137 >      constraintElems_.push_back(cp);
138      }
139 <
140 < }
141 <
121 < void Molecule::addConstraintElem(ConstraintElem* cp) {
122 <    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) {
123 <        constraintElems_.push_back(cp);
124 <    }
125 <
126 < }
127 <
128 < void Molecule::complete() {
139 >  }
140 >  
141 >  void Molecule::complete() {
142      
143      std::set<Atom*> rigidAtoms;
144 +    Atom* atom;
145 +    AtomIterator ai;
146      RigidBody* rb;
147 <    std::vector<RigidBody*>::iterator rbIter;
147 >    RigidBodyIterator rbIter;
148  
149 <    
149 >    // Get list of all the atoms that are part of rigid bodies
150 >
151      for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
152 <        rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
152 >      rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
153      }
154 +    
155 +    // add any atom that wasn't part of a rigid body to the list of integrableObjects
156  
139    Atom* atom;
140    AtomIterator ai;
157      for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
158 <  
159 <        if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
144 <            //if an atom does not belong to a rigid body, it is an integrable object
145 <            integrableObjects_.push_back(*ai);
146 <        }
147 <    }
158 >      
159 >      if (rigidAtoms.find(atom) == rigidAtoms.end()) {
160  
161 <    //find all free atoms (which do not belong to rigid bodies)  
162 <    //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_));
161 >        // If an atom does not belong to a rigid body, it is an
162 >        // integrable object
163  
164 <    //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
165 <    //    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
166 <    //    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
167 <    //
168 <    //    painCave.isFatal = 1;
169 <    //    simError();        
162 <    //}
164 >        integrableObjects_.push_back(atom);
165 >      }
166 >    }
167 >    
168 >    // then add the rigid bodies themselves to the integrableObjects
169 >
170      for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
171        integrableObjects_.push_back(rb);
172      }
166    //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
167 }
173  
174 < double Molecule::getMass() {
175 <    StuntDouble* sd;
171 <    std::vector<StuntDouble*>::iterator i;
172 <    double mass = 0.0;
174 >    // find the atoms that are fluctuating charges and add them to the
175 >    // fluctuatingCharges_ vector
176  
177 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
178 <        mass += sd->getMass();
177 >    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
178 >      if ( atom->isFluctuatingCharge() )
179 >        fluctuatingCharges_.push_back( atom );      
180      }
181  
182 <    return mass;
182 >  }
183  
184 < }
184 >  RealType Molecule::getMass() {
185 >    StuntDouble* sd;
186 >    std::vector<StuntDouble*>::iterator i;
187 >    RealType mass = 0.0;
188 >    
189 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
190 >           nextIntegrableObject(i)){
191 >      mass += sd->getMass();
192 >    }
193 >    
194 >    return mass;    
195 >  }
196  
197 < Vector3d Molecule::getCom() {
197 >  Vector3d Molecule::getCom() {
198      StuntDouble* sd;
199      std::vector<StuntDouble*>::iterator i;
200      Vector3d com;
201 <    double totalMass = 0;
202 <    double mass;
201 >    RealType totalMass = 0;
202 >    RealType mass;
203      
204 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
205 <        mass = sd->getMass();
206 <        totalMass += mass;
207 <        com += sd->getPos() * mass;    
204 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
205 >           nextIntegrableObject(i)){
206 >      mass = sd->getMass();
207 >      totalMass += mass;
208 >      com += sd->getPos() * mass;    
209      }
210 <
210 >    
211      com /= totalMass;
212  
213      return com;
214 < }
214 >  }
215  
216 < void Molecule::moveCom(const Vector3d& delta) {
216 >  void Molecule::moveCom(const Vector3d& delta) {
217      StuntDouble* sd;
218      std::vector<StuntDouble*>::iterator i;
219      
220 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
221 <        sd->setPos(sd->getPos() + delta);
222 <    }
220 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
221 >           nextIntegrableObject(i)){
222 >      sd->setPos(sd->getPos() + delta);
223 >    }    
224 >  }
225  
226 < }
209 <
210 < Vector3d Molecule::getComVel() {
226 >  Vector3d Molecule::getComVel() {
227      StuntDouble* sd;
228      std::vector<StuntDouble*>::iterator i;
229      Vector3d velCom;
230 <    double totalMass = 0;
231 <    double mass;
230 >    RealType totalMass = 0;
231 >    RealType mass;
232      
233 <    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
234 <        mass = sd->getMass();
235 <        totalMass += mass;
236 <        velCom += sd->getVel() * mass;    
233 >    for (sd = beginIntegrableObject(i); sd != NULL; sd =
234 >           nextIntegrableObject(i)){
235 >      mass = sd->getMass();
236 >      totalMass += mass;
237 >      velCom += sd->getVel() * mass;    
238      }
239 <
239 >    
240      velCom /= totalMass;
241 <
241 >    
242      return velCom;
243 < }
243 >  }
244  
245 < double Molecule::getPotential() {
245 >  RealType Molecule::getPotential() {
246  
247      Bond* bond;
248      Bend* bend;
249      Torsion* torsion;
250 +    Inversion* inversion;
251      Molecule::BondIterator bondIter;;
252      Molecule::BendIterator  bendIter;
253      Molecule::TorsionIterator  torsionIter;
254 +    Molecule::InversionIterator  inversionIter;
255  
256 <    double potential = 0.0;
256 >    RealType potential = 0.0;
257  
258      for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
259 <        potential += bond->getPotential();
259 >      potential += bond->getPotential();
260      }
261  
262      for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
263 <        potential += bend->getPotential();
263 >      potential += bend->getPotential();
264      }
265  
266 <    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
267 <        potential += torsion->getPotential();
266 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
267 >           nextTorsion(torsionIter)) {
268 >      potential += torsion->getPotential();
269      }
270 <
270 >    
271 >    for (inversion = beginInversion(inversionIter); torsion != NULL;
272 >         inversion =  nextInversion(inversionIter)) {
273 >      potential += inversion->getPotential();
274 >    }
275 >    
276      return potential;
277 +    
278 +  }
279 +  
280 +  void Molecule::addProperty(GenericData* genData) {
281 +    properties_.addProperty(genData);  
282 +  }
283  
284 < }
284 >  void Molecule::removeProperty(const std::string& propName) {
285 >    properties_.removeProperty(propName);  
286 >  }
287  
288 < std::ostream& operator <<(std::ostream& o, Molecule& mol) {
288 >  void Molecule::clearProperties() {
289 >    properties_.clearProperties();
290 >  }
291 >
292 >  std::vector<std::string> Molecule::getPropertyNames() {
293 >    return properties_.getPropertyNames();  
294 >  }
295 >      
296 >  std::vector<GenericData*> Molecule::getProperties() {
297 >    return properties_.getProperties();
298 >  }
299 >
300 >  GenericData* Molecule::getPropertyByName(const std::string& propName) {
301 >    return properties_.getPropertyByName(propName);
302 >  }
303 >
304 >  std::ostream& operator <<(std::ostream& o, Molecule& mol) {
305      o << std::endl;
306      o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
307      o << mol.getNAtoms() << " atoms" << std::endl;
308      o << mol.getNBonds() << " bonds" << std::endl;
309      o << mol.getNBends() << " bends" << std::endl;
310      o << mol.getNTorsions() << " torsions" << std::endl;
311 +    o << mol.getNInversions() << " inversions" << std::endl;
312      o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
313 <    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
314 <    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
315 <    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
313 >    o << mol.getNIntegrableObjects() << " integrable objects" << std::endl;
314 >    o << mol.getNCutoffGroups() << " cutoff groups" << std::endl;
315 >    o << mol.getNConstraintPairs() << " constraint pairs" << std::endl;
316 >    o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl;
317      return o;
318 < }
319 <
320 < }//end namespace oopse
318 >  }
319 >  
320 > }//end namespace OpenMD

Comparing:
trunk/src/primitives/Molecule.cpp (property svn:keywords), Revision 372 by tim, Mon Feb 21 15:28:56 2005 UTC vs.
branches/development/src/primitives/Molecule.cpp (property svn:keywords), Revision 1718 by gezelter, Thu May 24 01:29:59 2012 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines