ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/primitives/Molecule.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/primitives/Molecule.cpp (file contents):
Revision 1683, Thu Oct 28 22:34:02 2004 UTC vs.
Revision 1695 by tim, Mon Nov 1 22:52:57 2004 UTC

# Line 1 | Line 1
1 < #include <stdlib.h>
1 > /*
2 > * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
3 > *
4 > * Contact: oopse@oopse.org
5 > *
6 > * This program is free software; you can redistribute it and/or
7 > * modify it under the terms of the GNU Lesser General Public License
8 > * as published by the Free Software Foundation; either version 2.1
9 > * of the License, or (at your option) any later version.
10 > * All we ask is that proper credit is given for our work, which includes
11 > * - but is not limited to - adding the above copyright notice to the beginning
12 > * of your source code files, and to any copyright notice that you may distribute
13 > * with programs based on this work.
14 > *
15 > * This program is distributed in the hope that it will be useful,
16 > * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 > * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 > * GNU Lesser General Public License for more details.
19 > *
20 > * You should have received a copy of the GNU Lesser General Public License
21 > * along with this program; if not, write to the Free Software
22 > * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23 > *
24 > */
25  
26 + /**
27 + * @file Molecule.cpp
28 + * @author    tlin
29 + * @date  10/28/2004
30 + * @version 1.0
31 + */
32  
33 + #include <algorithm>
34 + #include <set>
35 +
36   #include "primitives/Molecule.hpp"
37 < #include "utils/simError.h"
37 > #include "utils/MemoryUtils.hpp"
38  
39 + namespace oopse {
40  
41 + Molecule::~Molecule() {
42  
43 < Molecule::Molecule( void ){
43 >    MemoryUtils::deleteVectorOfPointer(atoms_);
44 >    MemoryUtils::deleteVectorOfPointer(bonds_);
45 >    MemoryUtils::deleteVectorOfPointer(bends_);
46 >    MemoryUtils::deleteVectorOfPointer(torsions_);
47 >    MemoryUtils::deleteVectorOfPointer(rigidBodies_);
48 >    MemoryUtils::deleteVectorOfPointer(cutoffGroups_);
49  
50 <  myAtoms = NULL;
51 <  myBonds = NULL;
52 <  myBends = NULL;
14 <  myTorsions = NULL;
50 >    //integrableObjects_ don't own the objects
51 >    integrableObjects_.clear();
52 >    
53   }
54  
55 < Molecule::~Molecule( void ){
56 <  int i;
57 <  CutoffGroup* cg;
58 <  vector<CutoffGroup*>::iterator iter;
59 <  
22 <  if( myAtoms != NULL ){
23 <    for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
24 <    delete[] myAtoms;
25 <  }
55 > void Molecule::addAtom(Atom* atom) {
56 >    if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
57 >        atoms_.push_back(atom);
58 >    }
59 > }
60  
61 <  if( myBonds != NULL ){
62 <    for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
63 <    delete[] myBonds;
64 <  }
61 > void Molecule::addBond(Bond* bond) {
62 >    if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
63 >        bonds_.push_back(bond);
64 >    }
65 > }
66  
67 <  if( myBends != NULL ){
68 <    for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
69 <    delete[] myBends;
70 <  }
67 > void Molecule::addBend(Bend* bend) {
68 >    if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
69 >        bends_.push_back(bend);
70 >    }
71 > }
72  
73 <  if( myTorsions != NULL ){
74 <    for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
75 <    delete[] myTorsions;
76 <  }
73 > void Molecule::addTorsion(Torsion* torsion) {
74 >    if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
75 >        torsions_.push_back(torsion);
76 >    }
77 > }
78  
79 <  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
80 <    delete cg;
81 <  myCutoffGroups.clear();
82 <  
79 > void Molecule::addRigidBody(RigidBody *rb) {
80 >    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
81 >        rigidBodies_.push_back(rb);
82 >    }
83   }
84  
85 + void Molecule::addCutoffGroup(CutoffGroup* cp) {
86 +    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
87 +        cutoffGroups_.push_back(cp);
88 +    }
89  
90 < void Molecule::initialize( molInit &theInit ){
90 > }
91  
92 <  CutoffGroup* curCutoffGroup;
93 <  vector<CutoffGroup*>::iterator iterCutoff;
94 <  Atom* cutoffAtom;
95 <  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;
92 > void Molecule::complete() {
93 >    
94 >    std::set<Atom*> allAtoms;
95 >    allAtoms.insert(atoms_.begin(), atoms_.end());
96  
97 <  myAtoms = theInit.myAtoms;
98 <  myBonds = theInit.myBonds;
99 <  myBends = theInit.myBends;
68 <  myTorsions = theInit.myTorsions;
69 <  myRigidBodies = theInit.myRigidBodies;
97 >    std::set<Atom*> rigidAtoms;
98 >    RigidBody* rb;
99 >    std::vector<RigidBody*>::iterator rbIter;
100  
101 <  myIntegrableObjects = theInit.myIntegrableObjects;
101 >    
102 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
103 >        rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
104 >    }
105  
106 <  for (int i = 0; i < myRigidBodies.size(); i++)
107 <      myRigidBodies[i]->calcRefCoords();
106 >    //find all free atoms (which do not belong to rigid bodies)  
107 >    //performs the "difference" operation from set theory,  the output range contains a copy of every
108 >    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
109 >    //[rigidAtoms.begin(), rigidAtoms.end()).
110 >    std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
111 >                            std::back_inserter(integrableObjects_));
112  
113 <  myCutoffGroups = theInit.myCutoffGroups;
114 <  nCutoffGroups = myCutoffGroups.size();
113 >    if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
114 >        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
115 >    }
116  
117 +    integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
118   }
119  
120 < void Molecule::calcForces( void ){
121 <  
122 <  int i;
123 <  double com[3];
120 > Atom* Molecule::beginAtom(std::vector<Atom*>::iterator& i) {
121 >    i = atoms_.begin();
122 >    return (i == atoms_.end()) ? NULL : *i;
123 > }
124  
125 <  for(i=0; i<myRigidBodies.size(); i++) {
126 <    myRigidBodies[i]->updateAtoms();
127 <  }
125 > Atom* Molecule::nextAtom(std::vector<Atom*>::iterator& i) {
126 >    ++i;
127 >    return (i == atoms_.end()) ? NULL : *i;    
128 > }
129  
130 <  for(i=0; i<nBonds; i++){
131 <    myBonds[i]->calc_forces();
132 <  }
130 > Bond* Molecule::beginBond(std::vector<Bond*>::iterator& i) {
131 >    i = bonds_.begin();
132 >    return (i == bonds_.end()) ? NULL : *i;
133 > }
134  
135 <  for(i=0; i<nBends; i++){
136 <    myBends[i]->calc_forces();
137 <  }
135 > Bond* Molecule::nextBond(std::vector<Bond*>::iterator& i) {
136 >    ++i;
137 >    return (i == bonds_.end()) ? NULL : *i;    
138  
139 <  for(i=0; i<nTorsions; i++){
99 <    myTorsions[i]->calc_forces();
100 <  }
139 > }
140  
102 –  // Rigid Body forces and torques are done after the fortran force loop
141  
142 + Bend* Molecule::beginBend(std::vector<Bend*>::iterator& i) {
143 +    i = bends_.begin();
144 +    return (i == bends_.end()) ? NULL : *i;
145   }
146  
147 + Bend* Molecule::nextBend(std::vector<Bend*>::iterator& i) {
148 +    ++i;
149 +    return (i == bends_.end()) ? NULL : *i;    
150 + }
151  
152 < double Molecule::getPotential( void ){
153 <  
154 <  int i;
155 <  double myPot = 0.0;
152 > Torsion* Molecule::beginTorsion(std::vector<Torsion*>::iterator& i) {
153 >    i = torsions_.begin();
154 >    return (i == torsions_.end()) ? NULL : *i;
155 > }
156  
157 <  for(i=0; i<myRigidBodies.size(); i++) {
158 <    myRigidBodies[i]->updateAtoms();
159 <  }
160 <  
116 <  for(i=0; i<nBonds; i++){
117 <    myPot += myBonds[i]->get_potential();
118 <  }
157 > Torsion* Molecule::nextTorsion(std::vector<Torsion*>::iterator& i) {
158 >    ++i;
159 >    return (i == torsions_.end()) ? NULL : *i;    
160 > }    
161  
162 <  for(i=0; i<nBends; i++){
163 <    myPot += myBends[i]->get_potential();
164 <  }
162 > RigidBody* Molecule::beginRigidBody(std::vector<RigidBody*>::iterator& i) {
163 >    i = rigidBodies_.begin();
164 >    return (i == rigidBodies_.end()) ? NULL : *i;
165 > }
166  
167 <  for(i=0; i<nTorsions; i++){
168 <    myPot += myTorsions[i]->get_potential();
169 <  }
127 <
128 <  return myPot;
167 > RigidBody* Molecule::nextRigidBody(std::vector<RigidBody*>::iterator& i) {
168 >    ++i;
169 >    return (i == rigidBodies_.end()) ? NULL : *i;    
170   }
171  
172 < void Molecule::printMe( void ){
173 <  
174 <  int i;
175 <
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 <  }
172 > StuntDouble* Molecule::beginIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
173 >    i = integrableObjects_.begin();
174 >    return (i == integrableObjects_.end()) ? NULL : *i;
175 > }
176  
177 + StuntDouble* Molecule::nextIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
178 +    ++i;
179 +    return (i == integrableObjects_.end()) ? NULL : *i;    
180 + }    
181 +
182 + CutoffGroup* Molecule::beginCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
183 +    i = cutoffGroups_.begin();
184 +    return (i == cutoffGroups_.end()) ? NULL : *i;
185   }
186  
187 < void Molecule::moveCOM(double delta[3]){
188 <  double aPos[3];
189 <  int i, j;
187 > CutoffGroup* Molecule::nextCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {            
188 >    ++i;
189 >    return (i == cutoffGroups_.end()) ? NULL : *i;    
190 > }
191  
192 <  for(i=0; i<myIntegrableObjects.size(); i++) {
193 <    if(myIntegrableObjects[i] != NULL ) {
194 <      
195 <      myIntegrableObjects[i]->getPos( aPos );
196 <      
197 <      for (j=0; j< 3; j++)
198 <        aPos[j] += delta[j];
192 > void Molecule::calcForces() {
193 >    RigidBody* rb;
194 >    Bond* bond;
195 >    Bend* bend;
196 >    Torsion* torsion;
197 >    std::vector<RigidBody*>::iterator rbIter;
198 >    std::vector<Bond*>::iterator bondIter;;
199 >    std::vector<Bend*>::iterator bendIter;
200 >    std::vector<Torsion*>::iterator torsionIter;
201  
202 <      myIntegrableObjects[i]->setPos( aPos );
202 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
203 >        rb->updateAtoms();
204      }
163 –  }
205  
206 <  for(i=0; i<myRigidBodies.size(); i++) {
206 >    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
207 >        bond->calcForce();
208 >    }
209  
210 <      myRigidBodies[i]->getPos( aPos );
210 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
211 >        bend->calcForce();
212 >    }
213  
214 <      for (j=0; j< 3; j++)
215 <        aPos[j] += delta[j];
171 <      
172 <      myRigidBodies[i]->setPos( aPos );
214 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
215 >        torsion->calcForce();
216      }
217 +    
218   }
219  
220 < void Molecule::atoms2rigidBodies( void ) {
221 <  int i;
222 <  for (i = 0; i < myRigidBodies.size(); i++) {
223 <    myRigidBodies[i]->calcForcesAndTorques();  
224 <  }
225 < }
220 > double Molecule::getPotential() {
221 >    //RigidBody* rb;
222 >    Bond* bond;
223 >    Bend* bend;
224 >    Torsion* torsion;
225 >    //std::vector<RigidBody*> rbIter;
226 >    std::vector<Bond*>::iterator bondIter;;
227 >    std::vector<Bend*>::iterator bendIter;
228 >    std::vector<Torsion*>::iterator torsionIter;
229  
230 < void Molecule::getCOM( double COM[3] ) {
230 >    double potential = 0;
231 >    
232 >    //for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
233 >    //    rb->updateAtoms();
234 >    //}
235  
236 <  double mass, mtot;
237 <  double aPos[3];
238 <  int i, j;
236 >    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
237 >        potential += bond->getPotential();
238 >    }
239  
240 <  for (j=0; j<3; j++)
241 <    COM[j] = 0.0;
240 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
241 >        potential += bend->getPotential();
242 >    }
243  
244 <  mtot   = 0.0;
244 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
245 >        potential += torsion->getPotential();
246 >    }
247 >    
248 >    return potential;
249 > }
250  
251 <  for (i=0; i < myIntegrableObjects.size(); i++) {
252 <    if (myIntegrableObjects[i] != NULL) {
251 > Vector3d Molecule::getCom() {
252 >    StuntDouble* sd;
253 >    std::vector<StuntDouble*>::iterator i;
254 >    Vector3d com;
255 >    double totalMass = 0;
256 >    double mass;
257 >    
258 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
259 >        mass = sd->getMass();
260 >        totalMass += mass;
261 >        com += sd->getPos() * mass;    
262 >    }
263  
264 <      mass = myIntegrableObjects[i]->getMass();
198 <      mtot   += mass;
199 <      
200 <      myIntegrableObjects[i]->getPos( aPos );
264 >    com /= totalMass;
265  
266 <      for( j = 0; j < 3; j++)
267 <        COM[j] += aPos[j] * mass;
266 >    return com;
267 > }
268  
269 + void Molecule::moveCom(const Vector3d& delta) {
270 +    StuntDouble* sd;
271 +    std::vector<StuntDouble*>::iterator i;
272 +    
273 +    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
274 +        sd->setPos(sd->getPos() + delta);
275      }
206 –  }
276  
208 –  for (j = 0; j < 3; j++)
209 –    COM[j] /= mtot;
277   }
278  
279 < double Molecule::getCOMvel( double COMvel[3] ) {
280 <
281 <  double mass, mtot;
282 <  double aVel[3];
283 <  int i, j;
284 <
285 <
286 <  for (j=0; j<3; j++)
287 <    COMvel[j] = 0.0;
288 <
289 <  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 <
279 > Vector3d Molecule::getComVel() {
280 >    StuntDouble* sd;
281 >    std::vector<StuntDouble*>::iterator i;
282 >    Vector3d velCom;
283 >    double totalMass = 0;
284 >    double mass;
285 >    
286 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
287 >        mass = sd->getMass();
288 >        totalMass += mass;
289 >        velCom += sd->getVel() * mass;    
290      }
236 –  }
291  
292 <  for (j=0; j<3; j++)
239 <    COMvel[j] /= mtot;
240 <
241 <  return mtot;
292 >    velCom /= totalMass;
293  
294 +    return velCom;
295   }
296  
297 < double Molecule::getTotalMass()
298 < {
297 > std::ostream& operator <<(std::ostream& o, Molecule& mol) {
298 >    o << std::endl;
299 >    o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
300 >    o << mol.getNAtoms() << " atoms" << std::endl;
301 >    o << mol.getNBonds() << " bonds" << std::endl;
302 >    o << mol.getNBends() << " bends" << std::endl;
303 >    o << mol.getNTorsions() << " torsions" << std::endl;
304 >    o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
305 >    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
306 >    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
307  
308 <  double totalMass;
249 <  
250 <  totalMass = 0;
251 <  for(int i =0; i < myIntegrableObjects.size(); i++){
252 <    totalMass += myIntegrableObjects[i]->getMass();
253 <  }
254 <
255 <  return totalMass;
308 >    return o;
309   }
310 +
311 + }//end namespace oopse

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines