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

Comparing branches/new_design/OOPSE-4/src/primitives/Molecule.cpp (file contents):
Revision 1691, Thu Oct 28 22:34:02 2004 UTC vs.
Revision 1692 by tim, Mon Nov 1 20:15:58 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 "primitives/Molecule.hpp"
34 < #include "utils/simError.h"
34 > #include <algorithm>
35  
36 + namespace oopse {
37  
38 + Molecule::~Molecule() {
39  
40 < Molecule::Molecule( void ){
40 >    deleteVectorOfPointer(atoms_);
41 >    deleteVectorOfPointer(bonds_);
42 >    deleteVectorOfPointer(bends_);
43 >    deleteVectorOfPointer(torsions_);
44 >    deleteVectorOfPointer(rigidbodies_);
45 >    deleteVectorOfPointer(cutoffGroups_);
46  
47 <  myAtoms = NULL;
48 <  myBonds = NULL;
13 <  myBends = NULL;
14 <  myTorsions = NULL;
47 >    integrableObjects_.clear();
48 >    
49   }
50  
51 < Molecule::~Molecule( void ){
52 <  int i;
53 <  CutoffGroup* cg;
54 <  vector<CutoffGroup*>::iterator iter;
55 <  
22 <  if( myAtoms != NULL ){
23 <    for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
24 <    delete[] myAtoms;
25 <  }
51 > void Molecule::addAtom(Atom* atom) {
52 >    if (atoms_.find(atom) == atoms_.end()) {
53 >        atoms_.push_back(atom);
54 >    }
55 > }
56  
57 <  if( myBonds != NULL ){
58 <    for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
59 <    delete[] myBonds;
60 <  }
57 > void Molecule::addBond(Bond* bond) {
58 >    if (bonds_.find(bond) == bonds_.end()) {
59 >        bonds_.push_back(bond);
60 >    }
61 > }
62  
63 <  if( myBends != NULL ){
64 <    for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
65 <    delete[] myBends;
66 <  }
63 > void Molecule::addBend(Bend* bend) {
64 >    if (bends_.find(bend) == bends_.end()) {
65 >        bends_.push_back(bend);
66 >    }
67 > }
68  
69 <  if( myTorsions != NULL ){
70 <    for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
71 <    delete[] myTorsions;
72 <  }
69 > void Molecule::addTorsion(Torsion* torsion) {
70 >    if (torsions_.find(torsion) == torsions_.end()) {
71 >        torsions_.push_back(torsion);
72 >    }
73 > }
74  
75 <  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
76 <    delete cg;
77 <  myCutoffGroups.clear();
78 <  
75 > void Molecule::addRigidBody(RigidBody *rb) {
76 >    if (rigidBodies_.find(bond) == bonds_.end()) {
77 >        rigidBodies_.push_back(rb);
78 >    }
79   }
80  
81 + void Molecule::addCutoffGroup(CutoffGroup* cp) {
82 +    if (cutoffGroups_.find(bond) == bonds_.end()) {
83 +        cutoffGroups_.push_back(cp);
84 +    }
85  
86 < void Molecule::initialize( molInit &theInit ){
86 > }
87  
88 <  CutoffGroup* curCutoffGroup;
89 <  vector<CutoffGroup*>::iterator iterCutoff;
90 <  Atom* cutoffAtom;
91 <  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;
88 > void Molecule::complete() {
89 >    
90 >    std::set<Atom*> allAtoms;
91 >    allAtoms.insert(atoms_.begin(), atoms_.end());
92  
93 <  myAtoms = theInit.myAtoms;
94 <  myBonds = theInit.myBonds;
95 <  myBends = theInit.myBends;
68 <  myTorsions = theInit.myTorsions;
69 <  myRigidBodies = theInit.myRigidBodies;
93 >    std::set<Atom*> rigidAtoms;
94 >    RigidBody* rb;
95 >    std::vector<RigidBody*> rbIter;
96  
97 <  myIntegrableObjects = theInit.myIntegrableObjects;
97 >    
98 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
99 >        rigidAtoms.insert(rb->beginAtomIter(), rb->endAtomIter());
100 >    }
101  
102 <  for (int i = 0; i < myRigidBodies.size(); i++)
103 <      myRigidBodies[i]->calcRefCoords();
102 >    //find all free atoms (which do not belong to rigid bodies)  
103 >    //performs the "difference" operation from set theory,  the output range contains a copy of every
104 >    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
105 >    //[rigidAtoms.begin(), rigidAtoms.end()).
106 >    std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
107 >                            std::back_inserter(integrableObjects_));
108  
109 <  myCutoffGroups = theInit.myCutoffGroups;
110 <  nCutoffGroups = myCutoffGroups.size();
109 >    if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
110 >        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
111 >    }
112  
113 +    integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
114   }
115  
116 < void Molecule::calcForces( void ){
117 <  
118 <  int i;
119 <  double com[3];
116 > Atom* Molecule::beginAtom(std::vector<Atom*>::iterator& i) {
117 >    i = atoms_.begin();
118 >    return (i == atoms_.end()) ? NULL : *i;
119 > }
120  
121 <  for(i=0; i<myRigidBodies.size(); i++) {
122 <    myRigidBodies[i]->updateAtoms();
123 <  }
121 > Atom* Molecule::nextAtom(std::vector<Atom*>::iterator& i) {
122 >    ++i;
123 >    return (i == atoms_.end()) ? NULL : *i;    
124 > }
125  
126 <  for(i=0; i<nBonds; i++){
127 <    myBonds[i]->calc_forces();
128 <  }
126 > Bond* Molecule::beginBond(std::vector<Bond*>::iterator& i) {
127 >    i = bonds_.begin();
128 >    return (i == bonds_.end()) ? NULL : *i;
129 > }
130  
131 <  for(i=0; i<nBends; i++){
132 <    myBends[i]->calc_forces();
133 <  }
131 > Bond* Molecule::nextBond(std::vector<Bond*>::iterator& i) {
132 >    ++i;
133 >    return (i == bonds_.end()) ? NULL : *i;    
134  
135 <  for(i=0; i<nTorsions; i++){
99 <    myTorsions[i]->calc_forces();
100 <  }
135 > }
136  
102 –  // Rigid Body forces and torques are done after the fortran force loop
137  
138 + Bend* Molecule::beginBend(std::vector<Bend*>::iterator& i) {
139 +    i = bends_.begin();
140 +    return (i == bends_.end()) ? NULL : *i;
141   }
142  
143 + Bend* Molecule::nextBend(std::vector<Bend*>::iterator& i) {
144 +    ++i;
145 +    return (i == bends_.end()) ? NULL : *i;    
146 + }
147  
148 < double Molecule::getPotential( void ){
149 <  
150 <  int i;
151 <  double myPot = 0.0;
148 > Torsion* Molecule::beginTorsion(std::vector<Torsion*>::iterator& i) {
149 >    i = torsions_.begin();
150 >    return (i == torsions_.end()) ? NULL : *i;
151 > }
152  
153 <  for(i=0; i<myRigidBodies.size(); i++) {
154 <    myRigidBodies[i]->updateAtoms();
155 <  }
156 <  
116 <  for(i=0; i<nBonds; i++){
117 <    myPot += myBonds[i]->get_potential();
118 <  }
153 > Torsion* Molecule::nextTorsion(std::vector<Torsion*>::iterator& i) {
154 >    ++i;
155 >    return (i == torsions_.end()) ? NULL : *i;    
156 > }    
157  
158 <  for(i=0; i<nBends; i++){
159 <    myPot += myBends[i]->get_potential();
160 <  }
158 > RigidBody* Molecule::beginRigidBody(std::vector<RigidBody*>::iterator& i) {
159 >    i = rigidBodies_.begin();
160 >    return (i == rigidBodies_.end()) ? NULL : *i;
161 > }
162  
163 <  for(i=0; i<nTorsions; i++){
164 <    myPot += myTorsions[i]->get_potential();
165 <  }
127 <
128 <  return myPot;
163 > RigidBody* Molecule::nextRigidBody(std::vector<RigidBody*>::iterator& i) {
164 >    ++i;
165 >    return (i == rigidBodies_.end()) ? NULL : *i;    
166   }
167  
168 < void Molecule::printMe( void ){
169 <  
170 <  int i;
171 <
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 <  }
168 > StuntDouble* Molecule::beginIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
169 >    i = integrableObjects_.begin();
170 >    return (i == integrableObjects_.end()) ? NULL : *i;
171 > }
172  
173 + StuntDouble* Molecule::nextIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
174 +    ++i;
175 +    return (i == integrableObjects_.end()) ? NULL : *i;    
176 + }    
177 +
178 + CutoffGroup* Molecule::beginCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
179 +    i = cutoffGroups_.begin();
180 +    return (i == cutoffGroups_.end()) ? NULL : *i;
181   }
182  
183 < void Molecule::moveCOM(double delta[3]){
184 <  double aPos[3];
185 <  int i, j;
183 > CutoffGroup* Molecule::nextCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {            
184 >    ++i;
185 >    return (i == cutoffGroups_.end()) ? NULL : *i;    
186 > }
187  
188 <  for(i=0; i<myIntegrableObjects.size(); i++) {
189 <    if(myIntegrableObjects[i] != NULL ) {
190 <      
191 <      myIntegrableObjects[i]->getPos( aPos );
192 <      
193 <      for (j=0; j< 3; j++)
194 <        aPos[j] += delta[j];
188 > void Molecule::calcForces() {
189 >    RigidBody* rb;
190 >    Bond* bond;
191 >    Bend* bend;
192 >    Torsion* torsion;
193 >    std::vector<RigidBody*> rbIter;
194 >    std::vector<Bond*> bondIter;;
195 >    std::vector<Bend*> bendIter;
196 >    std::vector<Torsion*> torsionIter;
197  
198 <      myIntegrableObjects[i]->setPos( aPos );
198 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
199 >        rb->updateAtoms();
200      }
163 –  }
201  
202 <  for(i=0; i<myRigidBodies.size(); i++) {
202 >    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
203 >        bond->calcForce();
204 >    }
205  
206 <      myRigidBodies[i]->getPos( aPos );
206 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
207 >        bend->calcForce();
208 >    }
209  
210 <      for (j=0; j< 3; j++)
211 <        aPos[j] += delta[j];
171 <      
172 <      myRigidBodies[i]->setPos( aPos );
210 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
211 >        torsion->calcForce();
212      }
213 +    
214   }
215  
216 < void Molecule::atoms2rigidBodies( void ) {
217 <  int i;
218 <  for (i = 0; i < myRigidBodies.size(); i++) {
219 <    myRigidBodies[i]->calcForcesAndTorques();  
220 <  }
221 < }
216 > double Molecule::getPotential() {
217 >    //RigidBody* rb;
218 >    Bond* bond;
219 >    Bend* bend;
220 >    Torsion* torsion;
221 >    //std::vector<RigidBody*> rbIter;
222 >    std::vector<Bond*> bondIter;;
223 >    std::vector<Bend*> bendIter;
224 >    std::vector<Torsion*> torsionIter;
225  
226 < void Molecule::getCOM( double COM[3] ) {
226 >    double potential = 0;
227 >    
228 >    //for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
229 >    //    rb->updateAtoms();
230 >    //}
231  
232 <  double mass, mtot;
233 <  double aPos[3];
234 <  int i, j;
232 >    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
233 >        potential += bond->getPotential();
234 >    }
235  
236 <  for (j=0; j<3; j++)
237 <    COM[j] = 0.0;
236 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
237 >        potential += bend->getPotential();
238 >    }
239  
240 <  mtot   = 0.0;
240 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
241 >        potential += torsion->getPotential();
242 >    }
243 >    
244 >    return potential;
245 > }
246  
247 <  for (i=0; i < myIntegrableObjects.size(); i++) {
248 <    if (myIntegrableObjects[i] != NULL) {
247 > Vector3d Molecule::getCom() {
248 >    StuntDouble* sd;
249 >    std::vector<StuntDouble*>::iterator i;
250 >    Vector3d com;
251 >    double totalMass = 0;
252 >    double mass;
253 >    
254 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
255 >        mass = sd->getMass();
256 >        totalMass += mass;
257 >        com += sd->getPos() * mass;    
258 >    }
259  
260 <      mass = myIntegrableObjects[i]->getMass();
198 <      mtot   += mass;
199 <      
200 <      myIntegrableObjects[i]->getPos( aPos );
260 >    com /= totalMass;
261  
262 <      for( j = 0; j < 3; j++)
263 <        COM[j] += aPos[j] * mass;
262 >    return com;
263 > }
264  
265 + void Molecule::moveCom(const Vetor3d& delta) {
266 +    StuntDouble* sd;
267 +    std::vector<StuntDouble*>::iterator i;
268 +    
269 +    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
270 +        s->setPos(sd->getPos() + delta);
271      }
206 –  }
272  
208 –  for (j = 0; j < 3; j++)
209 –    COM[j] /= mtot;
273   }
274  
275 < double Molecule::getCOMvel( double COMvel[3] ) {
276 <
277 <  double mass, mtot;
278 <  double aVel[3];
279 <  int i, j;
280 <
281 <
282 <  for (j=0; j<3; j++)
283 <    COMvel[j] = 0.0;
284 <
285 <  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 <
275 > Vector3d Molecule::getComVel() {
276 >    StuntDouble* sd;
277 >    std::vector<StuntDouble*>::iterator i;
278 >    Vector3d velCom;
279 >    double totalMass = 0;
280 >    double mass;
281 >    
282 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
283 >        mass = sd->getMass();
284 >        totalMass += mass;
285 >        velCom += sd->getVel() * mass;    
286      }
236 –  }
287  
288 <  for (j=0; j<3; j++)
239 <    COMvel[j] /= mtot;
240 <
241 <  return mtot;
288 >    velCom /= totalMass;
289  
290 +    return velCom;
291   }
292  
293 < double Molecule::getTotalMass()
294 < {
293 > std::ostream& operator <<(std::ostream& o, const Molecule& mol) {
294 >    o << std::endl;
295 >    o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
296 >    o << mol.getNAtoms() << " atoms" << std::endl;
297 >    o << mol.getNBonds() << " bonds" << std::endl;
298 >    o << mol.getNBends() << " bends" << std::endl;
299 >    o << mol.getNTorsions() << " torsions" << std::endl;
300 >    o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
301 >    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
302 >    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
303  
304 <  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;
304 >    return o;
305   }
306 +
307 + }//end namespace oopse

Diff Legend

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