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

Diff Legend

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