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

Diff Legend

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