# | 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. 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 |
19 | > | * notice, this list of conditions and the following disclaimer. |
20 | > | * |
21 | > | * 3. Redistributions in binary form must reproduce the above copyright |
22 | > | * notice, this list of conditions and the following disclaimer in the |
23 | > | * documentation and/or other materials provided with the |
24 | > | * distribution. |
25 | > | * |
26 | > | * This software is provided "AS IS," without a warranty of any |
27 | > | * kind. All express or implied conditions, representations and |
28 | > | * warranties, including any implied warranty of merchantability, |
29 | > | * fitness for a particular purpose or non-infringement, are hereby |
30 | > | * excluded. The University of Notre Dame and its licensors shall not |
31 | > | * be liable for any damages suffered by licensee as a result of |
32 | > | * using, modifying or distributing the software or its |
33 | > | * derivatives. In no event will the University of Notre Dame or its |
34 | > | * licensors be liable for any lost revenue, profit or data, or for |
35 | > | * direct, indirect, special, consequential, incidental or punitive |
36 | > | * damages, however caused and regardless of the theory of liability, |
37 | > | * arising out of the use of or inability to use software, even if the |
38 | > | * University of Notre Dame has been advised of the possibility of |
39 | > | * such damages. |
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 "Molecule.hpp" |
53 | < | #include "simError.h" |
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 oopse { |
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(rigidBodies_); |
68 | > | MemoryUtils::deletePointers(cutoffGroups_); |
69 | > | MemoryUtils::deletePointers(constraintPairs_); |
70 | > | MemoryUtils::deletePointers(constraintElems_); |
71 | > | // integrableObjects_ don't own the objects |
72 | > | integrableObjects_.clear(); |
73 | > | |
74 | } | |
75 | < | |
76 | < | if( myBonds != NULL ){ |
77 | < | for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i]; |
78 | < | delete[] myBonds; |
75 | > | |
76 | > | void Molecule::addAtom(Atom* atom) { |
77 | > | if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) { |
78 | > | atoms_.push_back(atom); |
79 | > | } |
80 | } | |
81 | < | |
82 | < | if( myBends != NULL ){ |
83 | < | for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i]; |
84 | < | delete[] myBends; |
81 | > | |
82 | > | void Molecule::addBond(Bond* bond) { |
83 | > | if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) { |
84 | > | bonds_.push_back(bond); |
85 | > | } |
86 | } | |
87 | < | |
88 | < | if( myTorsions != NULL ){ |
89 | < | for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i]; |
90 | < | delete[] myTorsions; |
87 | > | |
88 | > | void Molecule::addBend(Bend* bend) { |
89 | > | if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) { |
90 | > | bends_.push_back(bend); |
91 | > | } |
92 | } | |
41 | – | |
42 | – | for(cg = beginCutoffGroup(iter); cg != NULL; cg = nextCutoffGroup(iter)) |
43 | – | delete cg; |
44 | – | myCutoffGroups.clear(); |
93 | ||
94 | < | } |
95 | < | |
96 | < | |
97 | < | void Molecule::initialize( molInit &theInit ){ |
98 | < | |
99 | < | CutoffGroup* curCutoffGroup; |
52 | < | vector<CutoffGroup*>::iterator iterCutoff; |
53 | < | Atom* cutoffAtom; |
54 | < | vector<Atom*>::iterator iterAtom; |
55 | < | int atomIndex; |
94 | > | void Molecule::addTorsion(Torsion* torsion) { |
95 | > | if (std::find(torsions_.begin(), torsions_.end(), torsion) == |
96 | > | torsions_.end()) { |
97 | > | torsions_.push_back(torsion); |
98 | > | } |
99 | > | } |
100 | ||
101 | < | nAtoms = theInit.nAtoms; |
102 | < | nMembers = nAtoms; |
103 | < | nBonds = theInit.nBonds; |
104 | < | nBends = theInit.nBends; |
105 | < | nTorsions = theInit.nTorsions; |
106 | < | 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 ){ |
101 | > | void Molecule::addRigidBody(RigidBody *rb) { |
102 | > | if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == |
103 | > | rigidBodies_.end()) { |
104 | > | rigidBodies_.push_back(rb); |
105 | > | } |
106 | > | } |
107 | ||
108 | < | int i; |
109 | < | double com[3]; |
110 | < | |
111 | < | for(i=0; i<myRigidBodies.size(); i++) { |
112 | < | myRigidBodies[i]->updateAtoms(); |
108 | > | void Molecule::addCutoffGroup(CutoffGroup* cp) { |
109 | > | if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == |
110 | > | cutoffGroups_.end()) { |
111 | > | cutoffGroups_.push_back(cp); |
112 | > | } |
113 | } | |
114 | < | |
115 | < | for(i=0; i<nBonds; i++){ |
116 | < | myBonds[i]->calc_forces(); |
114 | > | |
115 | > | void Molecule::addConstraintPair(ConstraintPair* cp) { |
116 | > | if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == |
117 | > | constraintPairs_.end()) { |
118 | > | constraintPairs_.push_back(cp); |
119 | > | } |
120 | } | |
121 | < | |
122 | < | for(i=0; i<nBends; i++){ |
123 | < | myBends[i]->calc_forces(); |
121 | > | |
122 | > | void Molecule::addConstraintElem(ConstraintElem* cp) { |
123 | > | if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == |
124 | > | constraintElems_.end()) { |
125 | > | constraintElems_.push_back(cp); |
126 | > | } |
127 | } | |
97 | – | |
98 | – | for(i=0; i<nTorsions; i++){ |
99 | – | myTorsions[i]->calc_forces(); |
100 | – | } |
101 | – | |
102 | – | // Rigid Body forces and torques are done after the fortran force loop |
103 | – | |
104 | – | } |
105 | – | |
106 | – | |
107 | – | double Molecule::getPotential( void ){ |
128 | ||
129 | < | int i; |
130 | < | double myPot = 0.0; |
129 | > | void Molecule::complete() { |
130 | > | |
131 | > | std::set<Atom*> rigidAtoms; |
132 | > | RigidBody* rb; |
133 | > | std::vector<RigidBody*>::iterator rbIter; |
134 | ||
135 | < | for(i=0; i<myRigidBodies.size(); i++) { |
136 | < | myRigidBodies[i]->updateAtoms(); |
137 | < | } |
138 | < | |
139 | < | for(i=0; i<nBonds; i++){ |
140 | < | myPot += myBonds[i]->get_potential(); |
141 | < | } |
135 | > | |
136 | > | for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { |
137 | > | rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); |
138 | > | } |
139 | > | |
140 | > | Atom* atom; |
141 | > | AtomIterator ai; |
142 | > | for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) { |
143 | > | |
144 | > | if (rigidAtoms.find(*ai) == rigidAtoms.end()) { |
145 | ||
146 | < | for(i=0; i<nBends; i++){ |
147 | < | myPot += myBends[i]->get_potential(); |
122 | < | } |
146 | > | // If an atom does not belong to a rigid body, it is an |
147 | > | // integrable object |
148 | ||
149 | < | for(i=0; i<nTorsions; i++){ |
150 | < | myPot += myTorsions[i]->get_potential(); |
151 | < | } |
149 | > | integrableObjects_.push_back(*ai); |
150 | > | } |
151 | > | } |
152 | > | |
153 | > | //find all free atoms (which do not belong to rigid bodies) |
154 | > | // performs the "difference" operation from set theory, the output |
155 | > | // range contains a copy of every element that is contained in |
156 | > | // [allAtoms.begin(), allAtoms.end()) and not contained in |
157 | > | // [rigidAtoms.begin(), rigidAtoms.end()). |
158 | > | //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(), |
159 | > | // std::back_inserter(integrableObjects_)); |
160 | ||
161 | < | return myPot; |
162 | < | } |
163 | < | |
164 | < | void Molecule::printMe( void ){ |
165 | < | |
166 | < | int i; |
167 | < | |
168 | < | for(i=0; i<nBonds; i++){ |
169 | < | myBonds[i]->printMe(); |
161 | > | //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) { |
162 | > | // //Some atoms in rigidAtoms are not in allAtoms, something must be wrong |
163 | > | // sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule"); |
164 | > | // |
165 | > | // painCave.isFatal = 1; |
166 | > | // simError(); |
167 | > | //} |
168 | > | for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { |
169 | > | integrableObjects_.push_back(rb); |
170 | > | } |
171 | > | //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end()); |
172 | } | |
173 | ||
174 | < | for(i=0; i<nBends; i++){ |
175 | < | myBends[i]->printMe(); |
176 | < | } |
177 | < | |
178 | < | for(i=0; i<nTorsions; i++){ |
179 | < | myTorsions[i]->printMe(); |
180 | < | } |
181 | < | |
147 | < | } |
148 | < | |
149 | < | void Molecule::moveCOM(double delta[3]){ |
150 | < | double aPos[3]; |
151 | < | int i, j; |
152 | < | |
153 | < | for(i=0; i<myIntegrableObjects.size(); i++) { |
154 | < | 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 ); |
174 | > | RealType Molecule::getMass() { |
175 | > | StuntDouble* sd; |
176 | > | std::vector<StuntDouble*>::iterator i; |
177 | > | RealType mass = 0.0; |
178 | > | |
179 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = |
180 | > | nextIntegrableObject(i)){ |
181 | > | mass += sd->getMass(); |
182 | } | |
183 | + | |
184 | + | return mass; |
185 | } | |
186 | ||
187 | < | for(i=0; i<myRigidBodies.size(); i++) { |
188 | < | |
189 | < | myRigidBodies[i]->getPos( aPos ); |
190 | < | |
191 | < | for (j=0; j< 3; j++) |
192 | < | aPos[j] += delta[j]; |
193 | < | |
194 | < | myRigidBodies[i]->setPos( aPos ); |
187 | > | Vector3d Molecule::getCom() { |
188 | > | StuntDouble* sd; |
189 | > | std::vector<StuntDouble*>::iterator i; |
190 | > | Vector3d com; |
191 | > | RealType totalMass = 0; |
192 | > | RealType mass; |
193 | > | |
194 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = |
195 | > | nextIntegrableObject(i)){ |
196 | > | mass = sd->getMass(); |
197 | > | totalMass += mass; |
198 | > | com += sd->getPos() * mass; |
199 | } | |
200 | < | } |
200 | > | |
201 | > | com /= totalMass; |
202 | ||
203 | < | void Molecule::atoms2rigidBodies( void ) { |
177 | < | int i; |
178 | < | for (i = 0; i < myRigidBodies.size(); i++) { |
179 | < | myRigidBodies[i]->calcForcesAndTorques(); |
203 | > | return com; |
204 | } | |
181 | – | } |
205 | ||
206 | < | void Molecule::getCOM( double COM[3] ) { |
206 | > | void Molecule::moveCom(const Vector3d& delta) { |
207 | > | StuntDouble* sd; |
208 | > | std::vector<StuntDouble*>::iterator i; |
209 | > | |
210 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = |
211 | > | nextIntegrableObject(i)){ |
212 | > | sd->setPos(sd->getPos() + delta); |
213 | > | } |
214 | > | } |
215 | ||
216 | < | double mass, mtot; |
217 | < | double aPos[3]; |
218 | < | int i, j; |
219 | < | |
220 | < | for (j=0; j<3; j++) |
221 | < | COM[j] = 0.0; |
222 | < | |
223 | < | mtot = 0.0; |
224 | < | |
225 | < | for (i=0; i < myIntegrableObjects.size(); i++) { |
226 | < | if (myIntegrableObjects[i] != NULL) { |
227 | < | |
197 | < | mass = myIntegrableObjects[i]->getMass(); |
198 | < | mtot += mass; |
199 | < | |
200 | < | myIntegrableObjects[i]->getPos( aPos ); |
201 | < | |
202 | < | for( j = 0; j < 3; j++) |
203 | < | COM[j] += aPos[j] * mass; |
204 | < | |
216 | > | Vector3d Molecule::getComVel() { |
217 | > | StuntDouble* sd; |
218 | > | std::vector<StuntDouble*>::iterator i; |
219 | > | Vector3d velCom; |
220 | > | RealType totalMass = 0; |
221 | > | RealType mass; |
222 | > | |
223 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = |
224 | > | nextIntegrableObject(i)){ |
225 | > | mass = sd->getMass(); |
226 | > | totalMass += mass; |
227 | > | velCom += sd->getVel() * mass; |
228 | } | |
229 | + | |
230 | + | velCom /= totalMass; |
231 | + | |
232 | + | return velCom; |
233 | } | |
234 | ||
235 | < | for (j = 0; j < 3; j++) |
209 | < | COM[j] /= mtot; |
210 | < | } |
235 | > | RealType Molecule::getPotential() { |
236 | ||
237 | < | double Molecule::getCOMvel( double COMvel[3] ) { |
237 | > | Bond* bond; |
238 | > | Bend* bend; |
239 | > | Torsion* torsion; |
240 | > | Molecule::BondIterator bondIter;; |
241 | > | Molecule::BendIterator bendIter; |
242 | > | Molecule::TorsionIterator torsionIter; |
243 | ||
244 | < | double mass, mtot; |
215 | < | double aVel[3]; |
216 | < | int i, j; |
244 | > | RealType potential = 0.0; |
245 | ||
246 | + | for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) { |
247 | + | potential += bond->getPotential(); |
248 | + | } |
249 | ||
250 | < | for (j=0; j<3; j++) |
251 | < | COMvel[j] = 0.0; |
250 | > | for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) { |
251 | > | potential += bend->getPotential(); |
252 | > | } |
253 | ||
254 | < | mtot = 0.0; |
255 | < | |
256 | < | 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 | < | |
254 | > | for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = |
255 | > | nextTorsion(torsionIter)) { |
256 | > | potential += torsion->getPotential(); |
257 | } | |
258 | + | |
259 | + | return potential; |
260 | + | |
261 | } | |
237 | – | |
238 | – | for (j=0; j<3; j++) |
239 | – | COMvel[j] /= mtot; |
240 | – | |
241 | – | return mtot; |
242 | – | |
243 | – | } |
244 | – | |
245 | – | double Molecule::getTotalMass() |
246 | – | { |
247 | – | |
248 | – | double totalMass; |
262 | ||
263 | < | totalMass = 0; |
264 | < | for(int i =0; i < myIntegrableObjects.size(); i++){ |
265 | < | totalMass += myIntegrableObjects[i]->getMass(); |
263 | > | std::ostream& operator <<(std::ostream& o, Molecule& mol) { |
264 | > | o << std::endl; |
265 | > | o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl; |
266 | > | o << mol.getNAtoms() << " atoms" << std::endl; |
267 | > | o << mol.getNBonds() << " bonds" << std::endl; |
268 | > | o << mol.getNBends() << " bends" << std::endl; |
269 | > | o << mol.getNTorsions() << " torsions" << std::endl; |
270 | > | o << mol.getNRigidBodies() << " rigid bodies" << std::endl; |
271 | > | o << mol.getNIntegrableObjects() << "integrable objects" << std::endl; |
272 | > | o << mol.getNCutoffGroups() << "cutoff groups" << std::endl; |
273 | > | o << mol.getNConstraintPairs() << "constraint pairs" << std::endl; |
274 | > | return o; |
275 | } | |
276 | < | |
277 | < | return totalMass; |
256 | < | } |
276 | > | |
277 | > | }//end namespace oopse |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |