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 483 by gezelter, Wed Apr 9 04:06:43 2003 UTC vs.
Revision 1158 by tim, Tue May 11 21:14:26 2004 UTC

# Line 1 | Line 1
1 < #include <cstdlib>
1 > #include <stdlib.h>
2  
3  
4   #include "Molecule.hpp"
# Line 12 | Line 12 | Molecule::Molecule( void ){
12    myBonds = NULL;
13    myBends = NULL;
14    myTorsions = NULL;
15 –
15   }
16  
18 –
19 –
17   Molecule::~Molecule( void ){
18    int i;
19    
# Line 40 | Line 37 | Molecule::~Molecule( void ){
37      delete[] myTorsions;
38    }
39  
43 –  if( myExcludes != NULL ){
44 –    for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
45 –    delete[] myExcludes;
46 –  }
40   }
41  
42  
43   void Molecule::initialize( molInit &theInit ){
44  
45 +  CutoffGroup* curCutoffGroup;
46 +  vector<CutoffGroup*>::iterator iterCutoff;
47 +  Atom* cutoffAtom;
48 +  vector<Atom*>::iterator iterAtom;
49 +  int atomIndex;
50 +  
51    nAtoms = theInit.nAtoms;
52    nMembers = nAtoms;
53    nBonds = theInit.nBonds;
54    nBends = theInit.nBends;
55    nTorsions = theInit.nTorsions;
56 <  nExcludes = theInit.nExcludes;
56 >  nRigidBodies = theInit.nRigidBodies;
57    nOriented = theInit.nOriented;
58 +  nCutoffGroups = theInit.nCutoffGroups;
59  
60    myAtoms = theInit.myAtoms;
61    myBonds = theInit.myBonds;
62    myBends = theInit.myBends;
63    myTorsions = theInit.myTorsions;
64 <  myExcludes = theInit.myExcludes;
65 <    
64 >  myRigidBodies = theInit.myRigidBodies;
65 >
66 >  myIntegrableObjects = theInit.myIntegrableObjects;
67 >
68 >  for (int i = 0; i < myRigidBodies.size(); i++)
69 >      myRigidBodies[i]->calcRefCoords();
70 >
71 >  myCutoffGroups = theInit.myCutoffGroups;
72 >
73 >  //creat a global index set of atoms which belong to cutoff group.
74 >  //it will be use to speed up the query whether an atom belongs to cutoff group or not
75 >  cutoffAtomSet.clear();
76 >
77 >  for(curCutoffGroup = beginCutoffGroup(iterCutoff); curCutoffGroup != NULL;
78 >                                                  curCutoffGroup = nextCutoffGroup(iterCutoff)){
79 >
80 >      for(cutoffAtom = curCutoffGroup->beginAtom(iterAtom); cutoffAtom != NULL;
81 >                                           cutoffAtom = curCutoffGroup->nextAtom(iterAtom)){
82 > #ifdef IS_MPI
83 >        atomIndex = cutoffAtom->getGlobalIndex();
84 > #else
85 >        atomIndex = cutoffAtom->getIndex();
86 > #endif
87 >        cutoffAtomSet.insert(atomIndex);
88 >     }// end for(cutoffAtom)    
89 >  }//end for(curCutoffGroup)
90 >  
91   }
92  
93   void Molecule::calcForces( void ){
94    
95    int i;
96 +  double com[3];
97  
98 +  for(i=0; i<myRigidBodies.size(); i++) {
99 +    myRigidBodies[i]->updateAtoms();
100 +  }
101 +
102 +  //calculate the center of mass of the molecule
103 +  //getCOM(com);  
104 +  //for(int i = 0; i < nAtoms; i ++)
105 +  //  myAtoms[i]->setRc(com);  
106 +  
107 +
108    for(i=0; i<nBonds; i++){
109      myBonds[i]->calc_forces();
110    }
# Line 80 | Line 116 | void Molecule::calcForces( void ){
116    for(i=0; i<nTorsions; i++){
117      myTorsions[i]->calc_forces();
118    }
119 +
120 +  // Rigid Body forces and torques are done after the fortran force loop
121 +
122   }
123  
124  
# Line 87 | Line 126 | double Molecule::getPotential( void ){
126    
127    int i;
128    double myPot = 0.0;
129 +
130 +  for(i=0; i<myRigidBodies.size(); i++) {
131 +    myRigidBodies[i]->updateAtoms();
132 +  }
133    
134    for(i=0; i<nBonds; i++){
135      myPot += myBonds[i]->get_potential();
# Line 118 | Line 161 | void Molecule::printMe( void ){
161    for(i=0; i<nTorsions; i++){
162      myTorsions[i]->printMe();
163    }
164 +
165   }
166  
167   void Molecule::moveCOM(double delta[3]){
168 <  double x, y, z;
169 <  int i;
168 >  double aPos[3];
169 >  int i, j;
170  
171 <  for(i=0; i<nAtoms; i++) {
172 <    if(myAtoms[i] != NULL ) {
171 >  for(i=0; i<myIntegrableObjects.size(); i++) {
172 >    if(myIntegrableObjects[i] != NULL ) {
173 >      
174 >      myIntegrableObjects[i]->getPos( aPos );
175 >      
176 >      for (j=0; j< 3; j++)
177 >        aPos[j] += delta[j];
178  
179 <      x = myAtoms[i]->getX() + delta[0];
180 <      y = myAtoms[i]->getY() + delta[1];
181 <      z = myAtoms[i]->getZ() + delta[2];
179 >      myIntegrableObjects[i]->setPos( aPos );
180 >    }
181 >  }
182  
183 <      myAtoms[i]->setX(x);
184 <      myAtoms[i]->setY(y);
185 <      myAtoms[i]->setZ(z);
183 >  for(i=0; i<myRigidBodies.size(); i++) {
184 >
185 >      myRigidBodies[i]->getPos( aPos );
186 >
187 >      for (j=0; j< 3; j++)
188 >        aPos[j] += delta[j];
189 >      
190 >      myRigidBodies[i]->setPos( aPos );
191      }
192 + }
193 +
194 + void Molecule::atoms2rigidBodies( void ) {
195 +  int i;
196 +  for (i = 0; i < myRigidBodies.size(); i++) {
197 +    myRigidBodies[i]->calcForcesAndTorques();  
198    }
199   }
200  
201   void Molecule::getCOM( double COM[3] ) {
202  
203    double mass, mtot;
204 <  int i;
204 >  double aPos[3];
205 >  int i, j;
206  
207 <  COM[0] = 0.0;
208 <  COM[1] = 0.0;
209 <  COM[2] = 0.0;
207 >  for (j=0; j<3; j++)
208 >    COM[j] = 0.0;
209 >
210    mtot   = 0.0;
211  
212 <  for (i=0; i < nAtoms; i++) {
213 <    if (myAtoms[i] != NULL) {
212 >  for (i=0; i < myIntegrableObjects.size(); i++) {
213 >    if (myIntegrableObjects[i] != NULL) {
214  
215 <      mass = myAtoms[i]->getMass();
215 >      mass = myIntegrableObjects[i]->getMass();
216        mtot   += mass;
217 <      COM[0] += myAtoms[i]->getX() * mass;
218 <      COM[1] += myAtoms[i]->getY() * mass;
158 <      COM[2] += myAtoms[i]->getZ() * mass;
217 >      
218 >      myIntegrableObjects[i]->getPos( aPos );
219  
220 +      for( j = 0; j < 3; j++)
221 +        COM[j] += aPos[j] * mass;
222 +
223      }
224    }
225  
226 <  COM[0] /= mtot;
227 <  COM[1] /= mtot;
165 <  COM[2] /= mtot;
166 <
226 >  for (j = 0; j < 3; j++)
227 >    COM[j] /= mtot;
228   }
229  
230   double Molecule::getCOMvel( double COMvel[3] ) {
231  
232    double mass, mtot;
233 <  int i;
233 >  double aVel[3];
234 >  int i, j;
235  
236 <  COMvel[0] = 0.0;
237 <  COMvel[1] = 0.0;
238 <  COMvel[2] = 0.0;
236 >
237 >  for (j=0; j<3; j++)
238 >    COMvel[j] = 0.0;
239 >
240    mtot   = 0.0;
241  
242 <  for (i=0; i < nAtoms; i++) {
243 <    if (myAtoms[i] != NULL) {
244 <      
245 <      mass = myAtoms[i]->getMass();
242 >  for (i=0; i < myIntegrableObjects.size(); i++) {
243 >    if (myIntegrableObjects[i] != NULL) {
244 >
245 >      mass = myIntegrableObjects[i]->getMass();
246        mtot   += mass;
184 –      COMvel[0] += myAtoms[i]->get_vx() * mass;
185 –      COMvel[1] += myAtoms[i]->get_vy() * mass;
186 –      COMvel[2] += myAtoms[i]->get_vz() * mass;
247  
248 +      myIntegrableObjects[i]->getVel(aVel);
249 +
250 +      for (j=0; j<3; j++)
251 +        COMvel[j] += aVel[j]*mass;
252 +
253      }
254    }
255  
256 <  COMvel[0] /= mtot;
257 <  COMvel[1] /= mtot;
193 <  COMvel[2] /= mtot;
256 >  for (j=0; j<3; j++)
257 >    COMvel[j] /= mtot;
258  
259    return mtot;
260  
261   }
262  
263 < void Molecule::atomicRollCall(int* molMembership) {
264 <  int i, which;
263 > double Molecule::getTotalMass()
264 > {
265  
266 <  for (i=0; i < nAtoms; i++) {
267 <    if (myAtoms[i] != NULL) {
268 <      which = myAtoms[i]->getIndex();
269 <      molMembership[which] = myIndex;
270 <    }
266 >  double totalMass;
267 >  
268 >  totalMass = 0;
269 >  for(int i =0; i < myIntegrableObjects.size(); i++){
270 >    totalMass += myIntegrableObjects[i]->getMass();
271    }
272 +
273 +  return totalMass;
274   }

Diff Legend

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