ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/integrators/Velocitizer.cpp
(Generate patch)

Comparing trunk/src/integrators/Velocitizer.cpp (file contents):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 555 by chuckv, Mon May 30 14:01:52 2005 UTC

# Line 50 | Line 50
50   #include "math/ParallelRandNumGen.hpp"
51   #endif
52  
53 + /* Remove me after testing*/
54 + #include <cstdio>
55 + #include <iostream>
56 + /*End remove me*/
57 +
58   namespace oopse {
59    
60    Velocitizer::Velocitizer(SimInfo* info) : info_(info) {
# Line 92 | Line 97 | namespace oopse {
97      double av2;
98      double kebar;
99      
100 +    Globals * simParams = info_->getSimParams();
101 +    
102      SimInfo::MoleculeIterator i;
103      Molecule::IntegrableObjectIterator j;
104      Molecule * mol;
# Line 149 | Line 156 | namespace oopse {
156      
157      
158      removeComDrift();
159 +    // Remove angular drift if we are not using periodic boundary conditions.
160 +    if(simParams->getPBC()) removeAngularDrift();
161      
162    }
163    
# Line 176 | Line 185 | namespace oopse {
185      
186    }
187    
188 +  
189 +   void Velocitizer::removeAngularDrift() {
190 +      // Get the Center of Mass drift velocity.
191 +      
192 +      Vector3d vdrift;
193 +      Vector3d com;
194 +      
195 +      info_->getComAll(com,vdrift);
196 +        
197 +      Mat3x3d inertiaTensor;
198 +      Vector3d angularMomentum;
199 +      Vector3d omega;
200 +      
201 +      
202 +      
203 +      info_->getInertiaTensor(inertiaTensor,angularMomentum);
204 +      // We now need the inverse of the inertia tensor.
205 +      
206 +      std::cerr << "Angular Momentum before is "
207 +         << angularMomentum <<  std::endl;
208 +
209 +      
210 +      inertiaTensor.inverse();
211 +      
212 +      
213 +      omega = inertiaTensor*angularMomentum;
214 +      
215 +      SimInfo::MoleculeIterator i;
216 +      Molecule::IntegrableObjectIterator j;
217 +      Molecule * mol;
218 +      StuntDouble * integrableObject;
219 +      Vector3d tempComPos;
220 +      
221 +      //  Corrects for the center of mass angular drift.
222 +      // sums all the angular momentum and divides by total mass.
223 +      for( mol = info_->beginMolecule(i); mol != NULL;
224 +           mol = info_->nextMolecule(i) ) {
225 +         for( integrableObject = mol->beginIntegrableObject(j);
226 +              integrableObject != NULL;
227 +              integrableObject = mol->nextIntegrableObject(j) ) {
228 +            tempComPos = integrableObject->getPos()-com;
229 +            integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos));
230 +         }
231 +      }
232 +      
233 +      angularMomentum = info_->getAngularMomentum();
234 +      std::cerr << "Angular Momentum after is "
235 +         << angularMomentum <<  std::endl;
236 +
237 +      
238 +   }
239 +  
240 +  
241 +  
242 +  
243   }

Diff Legend

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