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 621 by chuckv, Thu Sep 22 00:04:56 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 +      std::cerr << "Inertia Tensor before is "
209 +                << inertiaTensor <<  std::endl;
210 +      */
211 +      
212 +      inertiaTensor =inertiaTensor.inverse();
213 +       std::cerr << "Inertia Tensor after inverse is "
214 +                << inertiaTensor <<  std::endl;
215 +      
216 +      omega = inertiaTensor*angularMomentum;
217 +      
218 +      SimInfo::MoleculeIterator i;
219 +      Molecule::IntegrableObjectIterator j;
220 +      Molecule * mol;
221 +      StuntDouble * integrableObject;
222 +      Vector3d tempComPos;
223 +      
224 +      //  Corrects for the center of mass angular drift.
225 +      // sums all the angular momentum and divides by total mass.
226 +      for( mol = info_->beginMolecule(i); mol != NULL;
227 +           mol = info_->nextMolecule(i) ) {
228 +         for( integrableObject = mol->beginIntegrableObject(j);
229 +              integrableObject != NULL;
230 +              integrableObject = mol->nextIntegrableObject(j) ) {
231 +            tempComPos = integrableObject->getPos()-com;
232 +            integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos));
233 +         }
234 +      }
235 +      
236 +      angularMomentum = info_->getAngularMomentum();
237 +      /*
238 +      std::cerr << "Angular Momentum after is "
239 +         << angularMomentum <<  std::endl;
240 +      */
241 +      
242 +   }
243 +  
244 +  
245 +  
246 +  
247   }

Diff Legend

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