# | 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 | + | namespace oopse { |
57 | + | Molecule::Molecule(int stampId, int globalIndex, const std::string& molName) |
58 | + | : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) { |
59 | ||
60 | + | } |
61 | ||
62 | < | Molecule::Molecule( void ){ |
62 | > | Molecule::~Molecule() { |
63 | ||
64 | < | myAtoms = NULL; |
65 | < | myBonds = NULL; |
66 | < | myBends = NULL; |
67 | < | myTorsions = NULL; |
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 |
73 | > | integrableObjects_.clear(); |
74 | > | |
75 | } | |
76 | ||
77 | < | Molecule::~Molecule( void ){ |
78 | < | int i; |
79 | < | CutoffGroup* cg; |
80 | < | vector<CutoffGroup*>::iterator iter; |
81 | < | |
22 | < | if( myAtoms != NULL ){ |
23 | < | for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i]; |
24 | < | delete[] myAtoms; |
25 | < | } |
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( myBonds != NULL ){ |
84 | < | for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i]; |
85 | < | delete[] myBonds; |
86 | < | } |
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( myBends != NULL ){ |
90 | < | for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i]; |
91 | < | delete[] myBends; |
92 | < | } |
89 | > | void Molecule::addBend(Bend* bend) { |
90 | > | if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) { |
91 | > | bends_.push_back(bend); |
92 | > | } |
93 | > | } |
94 | ||
95 | < | if( myTorsions != NULL ){ |
96 | < | for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i]; |
97 | < | delete[] myTorsions; |
98 | < | } |
95 | > | void Molecule::addTorsion(Torsion* torsion) { |
96 | > | if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) { |
97 | > | torsions_.push_back(torsion); |
98 | > | } |
99 | > | } |
100 | ||
101 | < | for(cg = beginCutoffGroup(iter); cg != NULL; cg = nextCutoffGroup(iter)) |
102 | < | delete cg; |
103 | < | myCutoffGroups.clear(); |
104 | < | |
101 | > | void Molecule::addRigidBody(RigidBody *rb) { |
102 | > | if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) { |
103 | > | rigidBodies_.push_back(rb); |
104 | > | } |
105 | } | |
106 | ||
107 | + | void Molecule::addCutoffGroup(CutoffGroup* cp) { |
108 | + | if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) { |
109 | + | cutoffGroups_.push_back(cp); |
110 | + | } |
111 | ||
112 | < | void Molecule::initialize( molInit &theInit ){ |
112 | > | } |
113 | ||
114 | < | CutoffGroup* curCutoffGroup; |
115 | < | vector<CutoffGroup*>::iterator iterCutoff; |
116 | < | Atom* cutoffAtom; |
117 | < | 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; |
114 | > | void Molecule::addConstraintPair(ConstraintPair* cp) { |
115 | > | if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) { |
116 | > | constraintPairs_.push_back(cp); |
117 | > | } |
118 | ||
119 | < | myAtoms = theInit.myAtoms; |
66 | < | myBonds = theInit.myBonds; |
67 | < | myBends = theInit.myBends; |
68 | < | myTorsions = theInit.myTorsions; |
69 | < | myRigidBodies = theInit.myRigidBodies; |
119 | > | } |
120 | ||
121 | < | myIntegrableObjects = theInit.myIntegrableObjects; |
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 | < | for (int i = 0; i < myRigidBodies.size(); i++) |
74 | < | myRigidBodies[i]->calcRefCoords(); |
126 | > | } |
127 | ||
128 | < | myCutoffGroups = theInit.myCutoffGroups; |
129 | < | nCutoffGroups = myCutoffGroups.size(); |
128 | > | void Molecule::complete() { |
129 | > | |
130 | > | std::set<Atom*> rigidAtoms; |
131 | > | RigidBody* rb; |
132 | > | std::vector<RigidBody*>::iterator rbIter; |
133 | ||
134 | < | } |
134 | > | |
135 | > | for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { |
136 | > | rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); |
137 | > | } |
138 | ||
139 | < | void Molecule::calcForces( void ){ |
140 | < | |
141 | < | int i; |
142 | < | double com[3]; |
139 | > | Atom* atom; |
140 | > | AtomIterator ai; |
141 | > | for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) { |
142 | > | |
143 | > | 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 | > | } |
148 | ||
149 | < | for(i=0; i<myRigidBodies.size(); i++) { |
150 | < | myRigidBodies[i]->updateAtoms(); |
151 | < | } |
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 | < | for(i=0; i<nBonds; i++){ |
157 | < | myBonds[i]->calc_forces(); |
158 | < | } |
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 | > | //} |
163 | > | for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { |
164 | > | integrableObjects_.push_back(rb); |
165 | > | } |
166 | > | //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end()); |
167 | > | } |
168 | ||
169 | < | for(i=0; i<nBends; i++){ |
170 | < | myBends[i]->calc_forces(); |
171 | < | } |
169 | > | double Molecule::getMass() { |
170 | > | StuntDouble* sd; |
171 | > | std::vector<StuntDouble*>::iterator i; |
172 | > | double mass = 0.0; |
173 | ||
174 | < | for(i=0; i<nTorsions; i++){ |
175 | < | myTorsions[i]->calc_forces(); |
176 | < | } |
174 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ |
175 | > | mass += sd->getMass(); |
176 | > | } |
177 | ||
178 | < | // Rigid Body forces and torques are done after the fortran force loop |
178 | > | return mass; |
179 | ||
180 | } | |
181 | ||
182 | + | Vector3d Molecule::getCom() { |
183 | + | StuntDouble* sd; |
184 | + | std::vector<StuntDouble*>::iterator i; |
185 | + | Vector3d com; |
186 | + | double totalMass = 0; |
187 | + | double mass; |
188 | + | |
189 | + | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ |
190 | + | mass = sd->getMass(); |
191 | + | totalMass += mass; |
192 | + | com += sd->getPos() * mass; |
193 | + | } |
194 | ||
195 | < | double Molecule::getPotential( void ){ |
108 | < | |
109 | < | int i; |
110 | < | double myPot = 0.0; |
195 | > | com /= totalMass; |
196 | ||
197 | < | for(i=0; i<myRigidBodies.size(); i++) { |
113 | < | myRigidBodies[i]->updateAtoms(); |
114 | < | } |
115 | < | |
116 | < | for(i=0; i<nBonds; i++){ |
117 | < | myPot += myBonds[i]->get_potential(); |
118 | < | } |
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; |
197 | > | return com; |
198 | } | |
199 | ||
200 | < | void Molecule::printMe( void ){ |
201 | < | |
202 | < | int i; |
200 | > | void Molecule::moveCom(const Vector3d& delta) { |
201 | > | StuntDouble* sd; |
202 | > | std::vector<StuntDouble*>::iterator i; |
203 | > | |
204 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ |
205 | > | sd->setPos(sd->getPos() + delta); |
206 | > | } |
207 | ||
135 | – | for(i=0; i<nBonds; i++){ |
136 | – | myBonds[i]->printMe(); |
137 | – | } |
138 | – | |
139 | – | for(i=0; i<nBends; i++){ |
140 | – | myBends[i]->printMe(); |
141 | – | } |
142 | – | |
143 | – | for(i=0; i<nTorsions; i++){ |
144 | – | myTorsions[i]->printMe(); |
145 | – | } |
146 | – | |
208 | } | |
209 | ||
210 | < | void Molecule::moveCOM(double delta[3]){ |
211 | < | double aPos[3]; |
212 | < | int i, j; |
213 | < | |
214 | < | for(i=0; i<myIntegrableObjects.size(); i++) { |
215 | < | if(myIntegrableObjects[i] != NULL ) { |
216 | < | |
217 | < | myIntegrableObjects[i]->getPos( aPos ); |
218 | < | |
219 | < | for (j=0; j< 3; j++) |
220 | < | aPos[j] += delta[j]; |
160 | < | |
161 | < | myIntegrableObjects[i]->setPos( aPos ); |
210 | > | Vector3d Molecule::getComVel() { |
211 | > | StuntDouble* sd; |
212 | > | std::vector<StuntDouble*>::iterator i; |
213 | > | Vector3d velCom; |
214 | > | double totalMass = 0; |
215 | > | double mass; |
216 | > | |
217 | > | for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ |
218 | > | mass = sd->getMass(); |
219 | > | totalMass += mass; |
220 | > | velCom += sd->getVel() * mass; |
221 | } | |
163 | – | } |
222 | ||
223 | < | for(i=0; i<myRigidBodies.size(); i++) { |
223 | > | velCom /= totalMass; |
224 | ||
225 | < | myRigidBodies[i]->getPos( aPos ); |
168 | < | |
169 | < | for (j=0; j< 3; j++) |
170 | < | aPos[j] += delta[j]; |
171 | < | |
172 | < | myRigidBodies[i]->setPos( aPos ); |
173 | < | } |
225 | > | return velCom; |
226 | } | |
227 | ||
228 | < | void Molecule::atoms2rigidBodies( void ) { |
177 | < | int i; |
178 | < | for (i = 0; i < myRigidBodies.size(); i++) { |
179 | < | myRigidBodies[i]->calcForcesAndTorques(); |
180 | < | } |
181 | < | } |
228 | > | double Molecule::getPotential() { |
229 | ||
230 | < | void Molecule::getCOM( double COM[3] ) { |
230 | > | Bond* bond; |
231 | > | Bend* bend; |
232 | > | Torsion* torsion; |
233 | > | Molecule::BondIterator bondIter;; |
234 | > | Molecule::BendIterator bendIter; |
235 | > | Molecule::TorsionIterator torsionIter; |
236 | ||
237 | < | double mass, mtot; |
186 | < | double aPos[3]; |
187 | < | int i, j; |
237 | > | double potential = 0.0; |
238 | ||
239 | < | for (j=0; j<3; j++) |
240 | < | COM[j] = 0.0; |
239 | > | for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) { |
240 | > | potential += bond->getPotential(); |
241 | > | } |
242 | ||
243 | < | mtot = 0.0; |
244 | < | |
194 | < | for (i=0; i < myIntegrableObjects.size(); i++) { |
195 | < | if (myIntegrableObjects[i] != NULL) { |
196 | < | |
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 | < | |
243 | > | for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) { |
244 | > | potential += bend->getPotential(); |
245 | } | |
206 | – | } |
246 | ||
247 | < | for (j = 0; j < 3; j++) |
248 | < | COM[j] /= mtot; |
210 | < | } |
211 | < | |
212 | < | double Molecule::getCOMvel( double COMvel[3] ) { |
213 | < | |
214 | < | double mass, mtot; |
215 | < | double aVel[3]; |
216 | < | int i, j; |
217 | < | |
218 | < | |
219 | < | 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 | < | |
247 | > | for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) { |
248 | > | potential += torsion->getPotential(); |
249 | } | |
236 | – | } |
250 | ||
251 | < | for (j=0; j<3; j++) |
239 | < | COMvel[j] /= mtot; |
240 | < | |
241 | < | return mtot; |
251 | > | return potential; |
252 | ||
253 | } | |
254 | ||
255 | < | double Molecule::getTotalMass() |
256 | < | { |
257 | < | |
258 | < | double totalMass; |
259 | < | |
260 | < | totalMass = 0; |
261 | < | for(int i =0; i < myIntegrableObjects.size(); i++){ |
262 | < | totalMass += myIntegrableObjects[i]->getMass(); |
263 | < | } |
264 | < | |
265 | < | return totalMass; |
255 | > | std::ostream& operator <<(std::ostream& o, Molecule& mol) { |
256 | > | o << std::endl; |
257 | > | o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl; |
258 | > | o << mol.getNAtoms() << " atoms" << std::endl; |
259 | > | o << mol.getNBonds() << " bonds" << std::endl; |
260 | > | o << mol.getNBends() << " bends" << std::endl; |
261 | > | o << mol.getNTorsions() << " torsions" << std::endl; |
262 | > | o << mol.getNRigidBodies() << " rigid bodies" << std::endl; |
263 | > | o << mol.getNIntegrableObjects() << "integrable objects" << std::endl; |
264 | > | o << mol.getNCutoffGroups() << "cutoff groups" << std::endl; |
265 | > | o << mol.getNConstraintPairs() << "constraint pairs" << std::endl; |
266 | > | return o; |
267 | } | |
268 | + | |
269 | + | }//end namespace oopse |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |