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 1157 by tim, Tue May 11 20:33:41 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    
# Line 39 | Line 37 | Molecule::~Molecule( void ){
37      delete[] myTorsions;
38    }
39  
42 –  if( myExcludes != NULL ){
43 –    for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
44 –    delete[] myExcludes;
45 –  }
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.myExcludses;
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->beginAtom(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 79 | 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  
125 < void Molecule::getPotential( void ){
125 > 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 101 | Line 145 | void Molecule::getPotential( void ){
145  
146    return myPot;
147   }
148 +
149 + void Molecule::printMe( void ){
150 +  
151 +  int i;
152 +
153 +  for(i=0; i<nBonds; i++){
154 +    myBonds[i]->printMe();
155 +  }
156 +
157 +  for(i=0; i<nBends; i++){
158 +    myBends[i]->printMe();
159 +  }
160 +
161 +  for(i=0; i<nTorsions; i++){
162 +    myTorsions[i]->printMe();
163 +  }
164 +
165 + }
166 +
167 + void Molecule::moveCOM(double delta[3]){
168 +  double aPos[3];
169 +  int i, j;
170 +
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 +      myIntegrableObjects[i]->setPos( aPos );
180 +    }
181 +  }
182 +
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 +  double aPos[3];
205 +  int i, j;
206 +
207 +  for (j=0; j<3; j++)
208 +    COM[j] = 0.0;
209 +
210 +  mtot   = 0.0;
211 +
212 +  for (i=0; i < myIntegrableObjects.size(); i++) {
213 +    if (myIntegrableObjects[i] != NULL) {
214 +
215 +      mass = myIntegrableObjects[i]->getMass();
216 +      mtot   += 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 +  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 +  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;
241 +
242 +  for (i=0; i < myIntegrableObjects.size(); i++) {
243 +    if (myIntegrableObjects[i] != NULL) {
244 +
245 +      mass = myIntegrableObjects[i]->getMass();
246 +      mtot   += 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 +  for (j=0; j<3; j++)
257 +    COMvel[j] /= mtot;
258 +
259 +  return mtot;
260 +
261 + }
262 +
263 + double Molecule::getTotalMass()
264 + {
265 +
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