ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Molecule.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Molecule.cpp (file contents):
Revision 427 by mmeineke, Thu Mar 27 20:48:37 2003 UTC vs.
Revision 1452 by tim, Mon Aug 23 15:11:36 2004 UTC

# Line 1 | Line 1
1 < #include <cstdlib>
1 > #include <stdlib.h>
2  
3  
4   #include "Molecule.hpp"
5 + #include "simError.h"
6  
7  
8  
# Line 11 | Line 12 | Molecule::Molecule( void ){
12    myBonds = NULL;
13    myBends = NULL;
14    myTorsions = NULL;
14 –
15   }
16  
17 –
18 –
17   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];
# Line 39 | Line 39 | Molecule::~Molecule( void ){
39      delete[] myTorsions;
40    }
41  
42 <  if( myExcludes != NULL ){
43 <    for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
44 <    delete[] myExcludes;
45 <  }
42 >  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
43 >    delete cg;
44 >  myCutoffGroups.clear();
45 >  
46   }
47  
48  
49   void Molecule::initialize( molInit &theInit ){
50  
51 +  CutoffGroup* curCutoffGroup;
52 +  vector<CutoffGroup*>::iterator iterCutoff;
53 +  Atom* cutoffAtom;
54 +  vector<Atom*>::iterator iterAtom;
55 +  int atomIndex;
56 +  GenericData* gdata;
57 +  ConsRbData* rbData;
58 +  RigidBody* oldRb;
59 +  
60    nAtoms = theInit.nAtoms;
61    nMembers = nAtoms;
62    nBonds = theInit.nBonds;
63    nBends = theInit.nBends;
64    nTorsions = theInit.nTorsions;
65 <  nExcludes = theInit.nExcludes;
65 >  nRigidBodies = theInit.nRigidBodies;
66    nOriented = theInit.nOriented;
67  
68    myAtoms = theInit.myAtoms;
69    myBonds = theInit.myBonds;
70    myBends = theInit.myBends;
71    myTorsions = theInit.myTorsions;
72 <  myExcludes = theInit.myExcludses;
72 >  myRigidBodies = theInit.myRigidBodies;
73 >
74 >  myIntegrableObjects = theInit.myIntegrableObjects;
75 >
76 >  for (int i = 0; i < myRigidBodies.size(); i++){
77 >      myRigidBodies[i]->calcRefCoords();
78 >    //just a quick hack
79      
80 +    gdata = myRigidBodies[i]->getProperty("OldState");
81 +    if(gdata != NULL){
82 +      rbData = dynamic_cast<ConsRbData*>(gdata);
83 +      if(rbData ==NULL)
84 +        cerr << "dynamic_cast to ConsRbData Error in Molecule::initialize()" << endl;
85 +      else{
86 +        oldRb = rbData->getData();
87 +        oldRb->calcRefCoords();
88 +      }
89 +    }//end if(gata != NULL)
90 +    
91 +  }//end for(int i = 0; i < myRigidBodies.size(); i++)
92 +
93 +  myCutoffGroups = theInit.myCutoffGroups;
94 +  nCutoffGroups = myCutoffGroups.size();
95 +
96 +  myConstraintPairs = theInit.myConstraintPairs;
97 +  
98   }
99  
100   void Molecule::calcForces( void ){
101    
102    int i;
103 +  double com[3];
104  
105 +  for(i=0; i<myRigidBodies.size(); i++) {
106 +    myRigidBodies[i]->updateAtoms();
107 +  }
108 +
109 +  //calculate the center of mass of the molecule
110 +  //getCOM(com);  
111 +  //for(int i = 0; i < nAtoms; i ++)
112 +  //  myAtoms[i]->setRc(com);  
113 +  
114 +
115    for(i=0; i<nBonds; i++){
116      myBonds[i]->calc_forces();
117    }
# Line 79 | Line 123 | void Molecule::calcForces( void ){
123    for(i=0; i<nTorsions; i++){
124      myTorsions[i]->calc_forces();
125    }
126 +
127 +  // Rigid Body forces and torques are done after the fortran force loop
128 +
129   }
130  
131  
132 < void Molecule::getPotential( void ){
132 > double Molecule::getPotential( void ){
133    
134    int i;
135    double myPot = 0.0;
136 +
137 +  for(i=0; i<myRigidBodies.size(); i++) {
138 +    myRigidBodies[i]->updateAtoms();
139 +  }
140    
141    for(i=0; i<nBonds; i++){
142      myPot += myBonds[i]->get_potential();
# Line 101 | Line 152 | void Molecule::getPotential( void ){
152  
153    return myPot;
154   }
155 +
156 + void Molecule::printMe( void ){
157 +  
158 +  int i;
159 +
160 +  for(i=0; i<nBonds; i++){
161 +    myBonds[i]->printMe();
162 +  }
163 +
164 +  for(i=0; i<nBends; i++){
165 +    myBends[i]->printMe();
166 +  }
167 +
168 +  for(i=0; i<nTorsions; i++){
169 +    myTorsions[i]->printMe();
170 +  }
171 +
172 + }
173 +
174 + void Molecule::moveCOM(double delta[3]){
175 +  double aPos[3];
176 +  int i, j;
177 +
178 +  for(i=0; i<myIntegrableObjects.size(); i++) {
179 +    if(myIntegrableObjects[i] != NULL ) {
180 +      
181 +      myIntegrableObjects[i]->getPos( aPos );
182 +      
183 +      for (j=0; j< 3; j++)
184 +        aPos[j] += delta[j];
185 +
186 +      myIntegrableObjects[i]->setPos( aPos );
187 +    }
188 +  }
189 +
190 +  for(i=0; i<myRigidBodies.size(); i++) {
191 +
192 +      myRigidBodies[i]->getPos( aPos );
193 +
194 +      for (j=0; j< 3; j++)
195 +        aPos[j] += delta[j];
196 +      
197 +      myRigidBodies[i]->setPos( aPos );
198 +    }
199 + }
200 +
201 + void Molecule::atoms2rigidBodies( void ) {
202 +  int i;
203 +  for (i = 0; i < myRigidBodies.size(); i++) {
204 +    myRigidBodies[i]->calcForcesAndTorques();  
205 +  }
206 + }
207 +
208 + void Molecule::getCOM( double COM[3] ) {
209 +
210 +  double mass, mtot;
211 +  double aPos[3];
212 +  int i, j;
213 +
214 +  for (j=0; j<3; j++)
215 +    COM[j] = 0.0;
216 +
217 +  mtot   = 0.0;
218 +
219 +  for (i=0; i < myIntegrableObjects.size(); i++) {
220 +    if (myIntegrableObjects[i] != NULL) {
221 +
222 +      mass = myIntegrableObjects[i]->getMass();
223 +      mtot   += mass;
224 +      
225 +      myIntegrableObjects[i]->getPos( aPos );
226 +
227 +      for( j = 0; j < 3; j++)
228 +        COM[j] += aPos[j] * mass;
229 +
230 +    }
231 +  }
232 +
233 +  for (j = 0; j < 3; j++)
234 +    COM[j] /= mtot;
235 + }
236 +
237 + double Molecule::getCOMvel( double COMvel[3] ) {
238 +
239 +  double mass, mtot;
240 +  double aVel[3];
241 +  int i, j;
242 +
243 +
244 +  for (j=0; j<3; j++)
245 +    COMvel[j] = 0.0;
246 +
247 +  mtot   = 0.0;
248 +
249 +  for (i=0; i < myIntegrableObjects.size(); i++) {
250 +    if (myIntegrableObjects[i] != NULL) {
251 +
252 +      mass = myIntegrableObjects[i]->getMass();
253 +      mtot   += mass;
254 +
255 +      myIntegrableObjects[i]->getVel(aVel);
256 +
257 +      for (j=0; j<3; j++)
258 +        COMvel[j] += aVel[j]*mass;
259 +
260 +    }
261 +  }
262 +
263 +  for (j=0; j<3; j++)
264 +    COMvel[j] /= mtot;
265 +
266 +  return mtot;
267 +
268 + }
269 +
270 + double Molecule::getTotalMass()
271 + {
272 +
273 +  double totalMass;
274 +  
275 +  totalMass = 0;
276 +  for(int i =0; i < myIntegrableObjects.size(); i++){
277 +    totalMass += myIntegrableObjects[i]->getMass();
278 +  }
279 +
280 +  return totalMass;
281 + }

Diff Legend

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