# | Line 1 | Line 1 | |
---|---|---|
1 | < | #include <stdlib.h> |
1 | > | /* |
2 | > | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 | > | * |
4 | > | * The University of Notre Dame grants you ("Licensee") a |
5 | > | * non-exclusive, royalty free, license to use, modify and |
6 | > | * redistribute this software in source and binary code form, provided |
7 | > | * that the following conditions are met: |
8 | > | * |
9 | > | * 1. Redistributions of source code must retain the above copyright |
10 | > | * notice, this list of conditions and the following disclaimer. |
11 | > | * |
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. |
16 | > | * |
17 | > | * This software is provided "AS IS," without a warranty of any |
18 | > | * kind. All express or implied conditions, representations and |
19 | > | * warranties, including any implied warranty of merchantability, |
20 | > | * fitness for a particular purpose or non-infringement, are hereby |
21 | > | * excluded. The University of Notre Dame and its licensors shall not |
22 | > | * be liable for any damages suffered by licensee as a result of |
23 | > | * using, modifying or distributing the software or its |
24 | > | * derivatives. In no event will the University of Notre Dame or its |
25 | > | * licensors be liable for any lost revenue, profit or data, or for |
26 | > | * direct, indirect, special, consequential, incidental or punitive |
27 | > | * damages, however caused and regardless of the theory of liability, |
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] Vardeman & Gezelter, in progress (2009). |
40 | > | */ |
41 | > | |
42 | > | /** |
43 | > | * @file Molecule.cpp |
44 | > | * @author tlin |
45 | > | * @date 10/28/2004 |
46 | > | * @version 1.0 |
47 | > | */ |
48 | ||
49 | + | #include <algorithm> |
50 | + | #include <set> |
51 | ||
52 | #include "primitives/Molecule.hpp" | |
53 | + | #include "utils/MemoryUtils.hpp" |
54 | #include "utils/simError.h" | |
55 | ||
56 | < | |
57 | < | |
58 | < | Molecule::Molecule( void ){ |
59 | < | |
11 | < | myAtoms = NULL; |
12 | < | myBonds = NULL; |
13 | < | myBends = NULL; |
14 | < | myTorsions = NULL; |
15 | < | } |
16 | < | |
17 | < | Molecule::~Molecule( void ){ |
18 | < | int i; |
19 | < | CutoffGroup* cg; |
20 | < | vector<CutoffGroup*>::iterator iter; |
56 | > | namespace OpenMD { |
57 | > | Molecule::Molecule(int stampId, int globalIndex, const std::string& molName) |
58 | > | : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) { |
59 | > | } |
60 | ||
61 | < | if( myAtoms != NULL ){ |
62 | < | for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i]; |
63 | < | delete[] myAtoms; |
61 | > | Molecule::~Molecule() { |
62 | > | |
63 | > | MemoryUtils::deletePointers(atoms_); |
64 | > | MemoryUtils::deletePointers(bonds_); |
65 | > | MemoryUtils::deletePointers(bends_); |
66 | > | MemoryUtils::deletePointers(torsions_); |
67 | > | MemoryUtils::deletePointers(inversions_); |
68 | > | MemoryUtils::deletePointers(rigidBodies_); |
69 | > | MemoryUtils::deletePointers(cutoffGroups_); |
70 | > | MemoryUtils::deletePointers(constraintPairs_); |
71 | > | MemoryUtils::deletePointers(constraintElems_); |
72 | > | // integrableObjects_ don't own the objects |
73 | > | integrableObjects_.clear(); |
74 | > | |
75 | } | |
76 | < | |
77 | < | if( myBonds != NULL ){ |
78 | < | for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i]; |
79 | < | delete[] myBonds; |
76 | > | |
77 | > | void Molecule::addAtom(Atom* atom) { |
78 | > | if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) { |
79 | > | atoms_.push_back(atom); |
80 | > | } |
81 | } | |
82 | < | |
83 | < | if( myBends != NULL ){ |
84 | < | for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i]; |
85 | < | delete[] myBends; |
82 | > | |
83 | > | void Molecule::addBond(Bond* bond) { |
84 | > | if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) { |
85 | > | bonds_.push_back(bond); |
86 | > | } |
87 | } | |
88 | < | |
89 | < | if( myTorsions != NULL ){ |
90 | < | for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i]; |
91 | < | delete[] myTorsions; |
88 | > | |
89 | > | void Molecule::addBend(Bend* bend) { |
90 | > | if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) { |
91 | > | bends_.push_back(bend); |
92 | > | } |
93 | } | |
41 | – | |
42 | – | for(cg = beginCutoffGroup(iter); cg != NULL; cg = nextCutoffGroup(iter)) |
43 | – | delete cg; |
44 | – | myCutoffGroups.clear(); |
94 | ||
95 | < | } |
96 | < | |
97 | < | |
98 | < | void Molecule::initialize( molInit &theInit ){ |
99 | < | |
51 | < | CutoffGroup* curCutoffGroup; |
52 | < | vector<CutoffGroup*>::iterator iterCutoff; |
53 | < | Atom* cutoffAtom; |
54 | < | vector<Atom*>::iterator iterAtom; |
55 | < | int atomIndex; |
56 | < | |
57 | < | nAtoms = theInit.nAtoms; |
58 | < | nMembers = nAtoms; |
59 | < | nBonds = theInit.nBonds; |
60 | < | nBends = theInit.nBends; |
61 | < | nTorsions = theInit.nTorsions; |
62 | < | nRigidBodies = theInit.nRigidBodies; |
63 | < | nOriented = theInit.nOriented; |
64 | < | |
65 | < | myAtoms = theInit.myAtoms; |
66 | < | myBonds = theInit.myBonds; |
67 | < | myBends = theInit.myBends; |
68 | < | myTorsions = theInit.myTorsions; |
69 | < | myRigidBodies = theInit.myRigidBodies; |
70 | < | |
71 | < | myIntegrableObjects = theInit.myIntegrableObjects; |
72 | < | |
73 | < | for (int i = 0; i < myRigidBodies.size(); i++) |
74 | < | myRigidBodies[i]->calcRefCoords(); |
75 | < | |
76 | < | myCutoffGroups = theInit.myCutoffGroups; |
77 | < | nCutoffGroups = myCutoffGroups.size(); |
78 | < | |
79 | < | } |
80 | < | |
81 | < | void Molecule::calcForces( void ){ |
82 | < | |
83 | < | int i; |
84 | < | double com[3]; |
85 | < | |
86 | < | for(i=0; i<myRigidBodies.size(); i++) { |
87 | < | myRigidBodies[i]->updateAtoms(); |
95 | > | void Molecule::addTorsion(Torsion* torsion) { |
96 | > | if (std::find(torsions_.begin(), torsions_.end(), torsion) == |
97 | > | torsions_.end()) { |
98 | > | torsions_.push_back(torsion); |
99 | > | } |
100 | } | |
101 | ||
102 | < | for(i=0; i<nBonds; i++){ |
103 | < | myBonds[i]->calc_forces(); |
102 | > | void Molecule::addInversion(Inversion* inversion) { |
103 | > | if (std::find(inversions_.begin(), inversions_.end(), inversion) == |
104 | > | inversions_.end()) { |
105 | > | inversions_.push_back(inversion); |
106 | > | } |
107 | } | |
108 | < | |
109 | < | for(i=0; i<nBends; i++){ |
110 | < | myBends[i]->calc_forces(); |
108 | > | |
109 | > | void Molecule::addRigidBody(RigidBody *rb) { |
110 | > | if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == |
111 | > | rigidBodies_.end()) { |
112 | > | rigidBodies_.push_back(rb); |
113 | > | } |
114 | } | |
115 | < | |
116 | < | for(i=0; i<nTorsions; i++){ |
117 | < | myTorsions[i]->calc_forces(); |
115 | > | |
116 | > | void Molecule::addCutoffGroup(CutoffGroup* cp) { |
117 | > | if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == |
118 | > | cutoffGroups_.end()) { |
119 | > | cutoffGroups_.push_back(cp); |
120 | > | } |
121 | } | |
101 | – | |
102 | – | // Rigid Body forces and torques are done after the fortran force loop |
103 | – | |
104 | – | } |
105 | – | |
106 | – | |
107 | – | double Molecule::getPotential( void ){ |
122 | ||
123 | < | int i; |
124 | < | double myPot = 0.0; |
125 | < | |
126 | < | for(i=0; i<myRigidBodies.size(); i++) { |
127 | < | myRigidBodies[i]->updateAtoms(); |
123 | > | void Molecule::addConstraintPair(ConstraintPair* cp) { |
124 | > | if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == |
125 | > | constraintPairs_.end()) { |
126 | > | constraintPairs_.push_back(cp); |
127 | > | } |
128 | } | |
129 | < | |
130 | < | for(i=0; i<nBonds; i++){ |
131 | < | myPot += myBonds[i]->get_potential(); |
129 | > | |
130 | > | void Molecule::addConstraintElem(ConstraintElem* cp) { |
131 | > | if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == |
132 | > | constraintElems_.end()) { |
133 | > | constraintElems_.push_back(cp); |
134 | > | } |
135 | } | |
119 | – | |
120 | – | for(i=0; i<nBends; i++){ |
121 | – | myPot += myBends[i]->get_potential(); |
122 | – | } |
123 | – | |
124 | – | for(i=0; i<nTorsions; i++){ |
125 | – | myPot += myTorsions[i]->get_potential(); |
126 | – | } |
127 | – | |
128 | – | return myPot; |
129 | – | } |
130 | – | |
131 | – | void Molecule::printMe( void ){ |
136 | ||
137 | < | int i; |
137 | > | void Molecule::complete() { |
138 | > | |
139 | > | std::set<Atom*> rigidAtoms; |
140 | > | RigidBody* rb; |
141 | > | std::vector<RigidBody*>::iterator rbIter; |
142 | ||
143 | < | for(i=0; i<nBonds; i++){ |
144 | < | myBonds[i]->printMe(); |
145 | < | } |
143 | > | |
144 | > | for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { |
145 | > | rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); |
146 | > | } |
147 | > | |
148 | > | Atom* atom; |
149 | > | AtomIterator ai; |
150 | > | for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) { |
151 | > | |
152 | > | if (rigidAtoms.find(*ai) == rigidAtoms.end()) { |
153 | ||
154 | < | for(i=0; i<nBends; i++){ |
155 | < | myBends[i]->printMe(); |
141 | < | } |
154 | > | // If an atom does not belong to a rigid body, it is an |
155 | > | // integrable object |
156 | ||
157 | < | for(i=0; i<nTorsions; i++){ |
158 | < | myTorsions[i]->printMe(); |
157 | > | integrableObjects_.push_back(*ai); |
158 | > | } |
159 | > | } |
160 | > | |
161 | > | //find all free atoms (which do not belong to rigid bodies) |
162 | > | // performs the "difference" operation from set theory, the output |
163 | > | // range contains a copy of every element that is contained in |
164 | > | // [allAtoms.begin(), allAtoms.end()) and not contained in |
165 | > | // [rigidAtoms.begin(), rigidAtoms.end()). |
166 | > | //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(), |
167 | > | // std::back_inserter(integrableObjects_)); |
168 | > | |
169 | > | //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) { |
170 | > | // //Some atoms in rigidAtoms are not in allAtoms, something must be wrong |
171 | > | // sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule"); |
172 | > | // |
173 | > | // painCave.isFatal = 1; |
174 | > | // simError(); |
175 | > | //} |
176 | > | for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { |
177 | > | integrableObjects_.push_back(rb); |
178 | > | } |
179 | > | //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end()); |
180 | } | |
181 | ||
182 | < | } |
183 | < | |
184 | < | void Molecule::moveCOM(double delta[3]){ |
185 | < | double aPos[3]; |
186 | < | int i, j; |
187 | < | |
188 | < | for(i=0; i<myIntegrableObjects.size(); i++) { |
189 | < | if(myIntegrableObjects[i] != NULL ) { |
155 | < | |
156 | < | myIntegrableObjects[i]->getPos( aPos ); |
157 | < | |
158 | < | for (j=0; j< 3; j++) |
159 | < | aPos[j] += delta[j]; |
160 | < | |
161 | < | myIntegrableObjects[i]->setPos( aPos ); |
182 | > | RealType Molecule::getMass() { |
183 | > | StuntDouble* sd; |
184 | > | std::vector<StuntDouble*>::iterator i; |
185 | > | RealType mass = 0.0; |
186 | > | |
187 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = |
188 | > | nextIntegrableObject(i)){ |
189 | > | mass += sd->getMass(); |
190 | } | |
191 | + | |
192 | + | return mass; |
193 | } | |
194 | ||
195 | < | for(i=0; i<myRigidBodies.size(); i++) { |
196 | < | |
197 | < | myRigidBodies[i]->getPos( aPos ); |
198 | < | |
199 | < | for (j=0; j< 3; j++) |
200 | < | aPos[j] += delta[j]; |
201 | < | |
202 | < | myRigidBodies[i]->setPos( aPos ); |
195 | > | Vector3d Molecule::getCom() { |
196 | > | StuntDouble* sd; |
197 | > | std::vector<StuntDouble*>::iterator i; |
198 | > | Vector3d com; |
199 | > | RealType totalMass = 0; |
200 | > | RealType mass; |
201 | > | |
202 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = |
203 | > | nextIntegrableObject(i)){ |
204 | > | mass = sd->getMass(); |
205 | > | totalMass += mass; |
206 | > | com += sd->getPos() * mass; |
207 | } | |
208 | < | } |
208 | > | |
209 | > | com /= totalMass; |
210 | ||
211 | < | void Molecule::atoms2rigidBodies( void ) { |
177 | < | int i; |
178 | < | for (i = 0; i < myRigidBodies.size(); i++) { |
179 | < | myRigidBodies[i]->calcForcesAndTorques(); |
211 | > | return com; |
212 | } | |
181 | – | } |
213 | ||
214 | < | void Molecule::getCOM( double COM[3] ) { |
214 | > | void Molecule::moveCom(const Vector3d& delta) { |
215 | > | StuntDouble* sd; |
216 | > | std::vector<StuntDouble*>::iterator i; |
217 | > | |
218 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = |
219 | > | nextIntegrableObject(i)){ |
220 | > | sd->setPos(sd->getPos() + delta); |
221 | > | } |
222 | > | } |
223 | ||
224 | < | double mass, mtot; |
225 | < | double aPos[3]; |
226 | < | int i, j; |
224 | > | Vector3d Molecule::getComVel() { |
225 | > | StuntDouble* sd; |
226 | > | std::vector<StuntDouble*>::iterator i; |
227 | > | Vector3d velCom; |
228 | > | RealType totalMass = 0; |
229 | > | RealType mass; |
230 | > | |
231 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = |
232 | > | nextIntegrableObject(i)){ |
233 | > | mass = sd->getMass(); |
234 | > | totalMass += mass; |
235 | > | velCom += sd->getVel() * mass; |
236 | > | } |
237 | > | |
238 | > | velCom /= totalMass; |
239 | > | |
240 | > | return velCom; |
241 | > | } |
242 | ||
243 | < | for (j=0; j<3; j++) |
190 | < | COM[j] = 0.0; |
243 | > | RealType Molecule::getPotential() { |
244 | ||
245 | < | mtot = 0.0; |
245 | > | Bond* bond; |
246 | > | Bend* bend; |
247 | > | Torsion* torsion; |
248 | > | Inversion* inversion; |
249 | > | Molecule::BondIterator bondIter;; |
250 | > | Molecule::BendIterator bendIter; |
251 | > | Molecule::TorsionIterator torsionIter; |
252 | > | Molecule::InversionIterator inversionIter; |
253 | ||
254 | < | for (i=0; i < myIntegrableObjects.size(); i++) { |
195 | < | if (myIntegrableObjects[i] != NULL) { |
254 | > | RealType potential = 0.0; |
255 | ||
256 | < | mass = myIntegrableObjects[i]->getMass(); |
257 | < | mtot += mass; |
258 | < | |
200 | < | myIntegrableObjects[i]->getPos( aPos ); |
256 | > | for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) { |
257 | > | potential += bond->getPotential(); |
258 | > | } |
259 | ||
260 | < | for( j = 0; j < 3; j++) |
261 | < | COM[j] += aPos[j] * mass; |
260 | > | for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) { |
261 | > | potential += bend->getPotential(); |
262 | > | } |
263 | ||
264 | + | for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = |
265 | + | nextTorsion(torsionIter)) { |
266 | + | potential += torsion->getPotential(); |
267 | } | |
268 | + | |
269 | + | for (inversion = beginInversion(inversionIter); torsion != NULL; |
270 | + | inversion = nextInversion(inversionIter)) { |
271 | + | potential += inversion->getPotential(); |
272 | + | } |
273 | + | |
274 | + | return potential; |
275 | + | |
276 | } | |
277 | + | |
278 | + | void Molecule::addProperty(GenericData* genData) { |
279 | + | properties_.addProperty(genData); |
280 | + | } |
281 | ||
282 | < | for (j = 0; j < 3; j++) |
283 | < | COM[j] /= mtot; |
284 | < | } |
282 | > | void Molecule::removeProperty(const std::string& propName) { |
283 | > | properties_.removeProperty(propName); |
284 | > | } |
285 | ||
286 | < | double Molecule::getCOMvel( double COMvel[3] ) { |
286 | > | void Molecule::clearProperties() { |
287 | > | properties_.clearProperties(); |
288 | > | } |
289 | ||
290 | < | double mass, mtot; |
291 | < | double aVel[3]; |
292 | < | int i, j; |
290 | > | std::vector<std::string> Molecule::getPropertyNames() { |
291 | > | return properties_.getPropertyNames(); |
292 | > | } |
293 | > | |
294 | > | std::vector<GenericData*> Molecule::getProperties() { |
295 | > | return properties_.getProperties(); |
296 | > | } |
297 | ||
298 | < | |
299 | < | for (j=0; j<3; j++) |
220 | < | COMvel[j] = 0.0; |
221 | < | |
222 | < | mtot = 0.0; |
223 | < | |
224 | < | for (i=0; i < myIntegrableObjects.size(); i++) { |
225 | < | if (myIntegrableObjects[i] != NULL) { |
226 | < | |
227 | < | mass = myIntegrableObjects[i]->getMass(); |
228 | < | mtot += mass; |
229 | < | |
230 | < | myIntegrableObjects[i]->getVel(aVel); |
231 | < | |
232 | < | for (j=0; j<3; j++) |
233 | < | COMvel[j] += aVel[j]*mass; |
234 | < | |
235 | < | } |
298 | > | GenericData* Molecule::getPropertyByName(const std::string& propName) { |
299 | > | return properties_.getPropertyByName(propName); |
300 | } | |
301 | ||
238 | – | for (j=0; j<3; j++) |
239 | – | COMvel[j] /= mtot; |
240 | – | |
241 | – | return mtot; |
302 | ||
243 | – | } |
303 | ||
245 | – | double Molecule::getTotalMass() |
246 | – | { |
304 | ||
305 | < | double totalMass; |
306 | < | |
307 | < | totalMass = 0; |
308 | < | for(int i =0; i < myIntegrableObjects.size(); i++){ |
309 | < | totalMass += myIntegrableObjects[i]->getMass(); |
305 | > | std::ostream& operator <<(std::ostream& o, Molecule& mol) { |
306 | > | o << std::endl; |
307 | > | o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl; |
308 | > | o << mol.getNAtoms() << " atoms" << std::endl; |
309 | > | o << mol.getNBonds() << " bonds" << std::endl; |
310 | > | o << mol.getNBends() << " bends" << std::endl; |
311 | > | o << mol.getNTorsions() << " torsions" << std::endl; |
312 | > | o << mol.getNInversions() << " inversions" << std::endl; |
313 | > | o << mol.getNRigidBodies() << " rigid bodies" << std::endl; |
314 | > | o << mol.getNIntegrableObjects() << "integrable objects" << std::endl; |
315 | > | o << mol.getNCutoffGroups() << "cutoff groups" << std::endl; |
316 | > | o << mol.getNConstraintPairs() << "constraint pairs" << std::endl; |
317 | > | return o; |
318 | } | |
319 | < | |
320 | < | return totalMass; |
256 | < | } |
319 | > | |
320 | > | }//end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |