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

Comparing trunk/OOPSE/libmdtools/Integrator.cpp (file contents):
Revision 726 by tim, Tue Aug 26 20:37:30 2003 UTC vs.
Revision 778 by mmeineke, Fri Sep 19 20:00:27 2003 UTC

# Line 153 | Line 153 | template<typename T> void Integrator<T>::integrate(voi
153    double sampleTime = info->sampleTime;
154    double statusTime = info->statusTime;
155    double thermalTime = info->thermalTime;
156 +  double resetTime = info->resetTime;
157  
158 +
159    double currSample;
160    double currThermal;
161    double currStatus;
162 <
162 >  double currReset;
163 >  
164    int calcPot, calcStress;
165    int isError;
166  
# Line 174 | Line 177 | template<typename T> void Integrator<T>::integrate(voi
177    // initialize the forces before the first step
178  
179    calcForce(1, 1);
180 <  // myFF->doForces(1,1);
178 <
180 >  
181    if (info->setTemp){
182      thermalize();
183    }
184  
183 –  calcPot = 0;
184 –  calcStress = 0;
185 –  currSample = sampleTime;
186 –  currThermal = thermalTime;
187 –  currStatus = statusTime;
188 –  
185    calcPot     = 0;
186    calcStress  = 0;
187    currSample  = sampleTime + info->getTime();
188    currThermal = thermalTime+ info->getTime();
189    currStatus  = statusTime + info->getTime();
190 +  currReset   = resetTime  + info->getTime();
191  
192    dumpOut->writeDump(info->getTime());
193    statOut->writeStat(info->getTime());
# Line 231 | Line 228 | template<typename T> void Integrator<T>::integrate(voi
228        currStatus += statusTime;
229      }
230  
231 +    if (info->resetIntegrator){
232 +      if (info->getTime() >= currReset){
233 +        this->resetIntegrator();
234 +        currReset += resetTime;
235 +      }
236 +    }
237 +
238   #ifdef IS_MPI
239      strcpy(checkPointMsg, "successfully took a time step.");
240      MPIcheckPoint();
# Line 250 | Line 254 | template<typename T> void Integrator<T>::integrateStep
254  
255    moveA();
256  
253 –  if (nConstrained){
254 –    constrainA();
255 –  }
257  
258  
259 +
260   #ifdef IS_MPI
261    strcpy(checkPointMsg, "Succesful moveA\n");
262    MPIcheckPoint();
# Line 275 | Line 277 | template<typename T> void Integrator<T>::integrateStep
277  
278    moveB();
279  
278 –  if (nConstrained){
279 –    constrainB();
280 –  }
280  
281 +
282   #ifdef IS_MPI
283    strcpy(checkPointMsg, "Succesful moveB\n");
284    MPIcheckPoint();
# Line 290 | Line 290 | template<typename T> void Integrator<T>::moveA(void){
290    int i, j;
291    DirectionalAtom* dAtom;
292    double Tb[3], ji[3];
293 –  double A[3][3], I[3][3];
294 –  double angle;
293    double vel[3], pos[3], frc[3];
294    double mass;
295  
# Line 327 | Line 325 | template<typename T> void Integrator<T>::moveA(void){
325        for (j = 0; j < 3; j++)
326          ji[j] += (dt2 * Tb[j]) * eConvert;
327  
328 <      // use the angular velocities to propagate the rotation matrix a
331 <      // full time step
328 >      this->rotationPropagation( dAtom, ji );
329  
330 <      dAtom->getA(A);
331 <      dAtom->getI(I);
332 <
336 <      // rotate about the x-axis      
337 <      angle = dt2 * ji[0] / I[0][0];
338 <      this->rotate(1, 2, angle, ji, A);
339 <
340 <      // rotate about the y-axis
341 <      angle = dt2 * ji[1] / I[1][1];
342 <      this->rotate(2, 0, angle, ji, A);
343 <
344 <      // rotate about the z-axis
345 <      angle = dt * ji[2] / I[2][2];
346 <      this->rotate(0, 1, angle, ji, A);
347 <
348 <      // rotate about the y-axis
349 <      angle = dt2 * ji[1] / I[1][1];
350 <      this->rotate(2, 0, angle, ji, A);
330 >      dAtom->setJ(ji);
331 >    }
332 >  }
333  
334 <      // rotate about the x-axis
335 <      angle = dt2 * ji[0] / I[0][0];
354 <      this->rotate(1, 2, angle, ji, A);
355 <
356 <
357 <      dAtom->setJ(ji);
358 <      dAtom->setA(A);
359 <    }
334 >  if (nConstrained){
335 >    constrainA();
336    }
337   }
338  
# Line 399 | Line 375 | template<typename T> void Integrator<T>::moveB(void){
375        dAtom->setJ(ji);
376      }
377    }
378 +
379 +  if (nConstrained){
380 +    constrainB();
381 +  }
382   }
383  
384   template<typename T> void Integrator<T>::preMove(void){
# Line 557 | Line 537 | template<typename T> void Integrator<T>::constrainA(){
537      painCave.isFatal = 1;
538      simError();
539    }
540 +
541   }
542  
543   template<typename T> void Integrator<T>::constrainB(void){
# Line 658 | Line 639 | template<typename T> void Integrator<T>::constrainB(vo
639      painCave.isFatal = 1;
640      simError();
641    }
642 + }
643 +
644 + template<typename T> void Integrator<T>::rotationPropagation
645 + ( DirectionalAtom* dAtom, double ji[3] ){
646 +
647 +  double angle;
648 +  double A[3][3], I[3][3];
649 +
650 +  // use the angular velocities to propagate the rotation matrix a
651 +  // full time step
652 +
653 +  dAtom->getA(A);
654 +  dAtom->getI(I);
655 +  
656 +  // rotate about the x-axis      
657 +  angle = dt2 * ji[0] / I[0][0];
658 +  this->rotate( 1, 2, angle, ji, A );
659 +  
660 +  // rotate about the y-axis
661 +  angle = dt2 * ji[1] / I[1][1];
662 +  this->rotate( 2, 0, angle, ji, A );
663 +  
664 +  // rotate about the z-axis
665 +  angle = dt * ji[2] / I[2][2];
666 +  this->rotate( 0, 1, angle, ji, A);
667 +  
668 +  // rotate about the y-axis
669 +  angle = dt2 * ji[1] / I[1][1];
670 +  this->rotate( 2, 0, angle, ji, A );
671 +  
672 +  // rotate about the x-axis
673 +  angle = dt2 * ji[0] / I[0][0];
674 +  this->rotate( 1, 2, angle, ji, A );
675 +  
676 +  dAtom->setA( A  );    
677   }
678  
679   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
# Line 750 | Line 766 | template<typename T> void Integrator<T>::thermalize(){
766   template<typename T> void Integrator<T>::thermalize(){
767    tStats->velocitize();
768   }
769 +
770 + template<typename T> double Integrator<T>::getConservedQuantity(void){
771 +  return tStats->getTotalE();
772 + }

Diff Legend

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