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 733 by tim, Wed Aug 27 19:23:29 2003 UTC vs.
Revision 782 by mmeineke, Tue Sep 23 20:34:31 2003 UTC

# Line 25 | Line 25 | template<typename T> Integrator<T>::Integrator(SimInfo
25    if (info->the_integrator != NULL){
26      delete info->the_integrator;
27    }
28 <  info->the_integrator = this;
29 <
28 >  
29    nAtoms = info->n_atoms;
30  
31    // check for constraints
# Line 153 | Line 152 | template<typename T> void Integrator<T>::integrate(voi
152    double sampleTime = info->sampleTime;
153    double statusTime = info->statusTime;
154    double thermalTime = info->thermalTime;
155 +  double resetTime = info->resetTime;
156  
157 +
158    double currSample;
159    double currThermal;
160    double currStatus;
161 <
161 >  double currReset;
162 >  
163    int calcPot, calcStress;
164    int isError;
165  
# Line 174 | Line 176 | template<typename T> void Integrator<T>::integrate(voi
176    // initialize the forces before the first step
177  
178    calcForce(1, 1);
179 +
180 +  if (nConstrained){
181 +    preMove();
182 +    constrainA();
183 +    calcForce(1, 1);    
184 +    constrainB();
185 +  }
186    
187    if (info->setTemp){
188      thermalize();
189    }
190  
182 –  calcPot = 0;
183 –  calcStress = 0;
184 –  currSample = sampleTime;
185 –  currThermal = thermalTime;
186 –  currStatus = statusTime;
187 –  
191    calcPot     = 0;
192    calcStress  = 0;
193    currSample  = sampleTime + info->getTime();
194    currThermal = thermalTime+ info->getTime();
195    currStatus  = statusTime + info->getTime();
196 +  currReset   = resetTime  + info->getTime();
197  
198    dumpOut->writeDump(info->getTime());
199    statOut->writeStat(info->getTime());
# Line 230 | Line 234 | template<typename T> void Integrator<T>::integrate(voi
234        currStatus += statusTime;
235      }
236  
237 +    if (info->resetIntegrator){
238 +      if (info->getTime() >= currReset){
239 +        this->resetIntegrator();
240 +        currReset += resetTime;
241 +      }
242 +    }
243 +
244   #ifdef IS_MPI
245      strcpy(checkPointMsg, "successfully took a time step.");
246      MPIcheckPoint();
# Line 249 | Line 260 | template<typename T> void Integrator<T>::integrateStep
260  
261    moveA();
262  
263 <  if (nConstrained){
253 <    constrainA();
254 <  }
263 >
264  
265  
266   #ifdef IS_MPI
# Line 274 | Line 283 | template<typename T> void Integrator<T>::integrateStep
283  
284    moveB();
285  
277 –  if (nConstrained){
278 –    constrainB();
279 –  }
286  
287 +
288   #ifdef IS_MPI
289    strcpy(checkPointMsg, "Succesful moveB\n");
290    MPIcheckPoint();
# Line 289 | Line 296 | template<typename T> void Integrator<T>::moveA(void){
296    int i, j;
297    DirectionalAtom* dAtom;
298    double Tb[3], ji[3];
292 –  double A[3][3], I[3][3];
293 –  double angle;
299    double vel[3], pos[3], frc[3];
300    double mass;
301  
# Line 326 | Line 331 | template<typename T> void Integrator<T>::moveA(void){
331        for (j = 0; j < 3; j++)
332          ji[j] += (dt2 * Tb[j]) * eConvert;
333  
334 <      // use the angular velocities to propagate the rotation matrix a
330 <      // full time step
334 >      this->rotationPropagation( dAtom, ji );
335  
336 <      dAtom->getA(A);
337 <      dAtom->getI(I);
336 >      dAtom->setJ(ji);
337 >    }
338 >  }
339  
340 <      // rotate about the x-axis      
341 <      angle = dt2 * ji[0] / I[0][0];
337 <      this->rotate(1, 2, angle, ji, A);
338 <
339 <      // rotate about the y-axis
340 <      angle = dt2 * ji[1] / I[1][1];
341 <      this->rotate(2, 0, angle, ji, A);
342 <
343 <      // rotate about the z-axis
344 <      angle = dt * ji[2] / I[2][2];
345 <      this->rotate(0, 1, angle, ji, A);
346 <
347 <      // rotate about the y-axis
348 <      angle = dt2 * ji[1] / I[1][1];
349 <      this->rotate(2, 0, angle, ji, A);
350 <
351 <      // rotate about the x-axis
352 <      angle = dt2 * ji[0] / I[0][0];
353 <      this->rotate(1, 2, angle, ji, A);
354 <
355 <
356 <      dAtom->setJ(ji);
357 <      dAtom->setA(A);
358 <    }
340 >  if (nConstrained){
341 >    constrainA();
342    }
343   }
344  
# Line 398 | Line 381 | template<typename T> void Integrator<T>::moveB(void){
381        dAtom->setJ(ji);
382      }
383    }
384 +
385 +  if (nConstrained){
386 +    constrainB();
387 +  }
388   }
389  
390   template<typename T> void Integrator<T>::preMove(void){
# Line 556 | Line 543 | template<typename T> void Integrator<T>::constrainA(){
543      painCave.isFatal = 1;
544      simError();
545    }
546 +
547   }
548  
549   template<typename T> void Integrator<T>::constrainB(void){
# Line 657 | Line 645 | template<typename T> void Integrator<T>::constrainB(vo
645      painCave.isFatal = 1;
646      simError();
647    }
648 + }
649 +
650 + template<typename T> void Integrator<T>::rotationPropagation
651 + ( DirectionalAtom* dAtom, double ji[3] ){
652 +
653 +  double angle;
654 +  double A[3][3], I[3][3];
655 +
656 +  // use the angular velocities to propagate the rotation matrix a
657 +  // full time step
658 +
659 +  dAtom->getA(A);
660 +  dAtom->getI(I);
661 +  
662 +  // rotate about the x-axis      
663 +  angle = dt2 * ji[0] / I[0][0];
664 +  this->rotate( 1, 2, angle, ji, A );
665 +  
666 +  // rotate about the y-axis
667 +  angle = dt2 * ji[1] / I[1][1];
668 +  this->rotate( 2, 0, angle, ji, A );
669 +  
670 +  // rotate about the z-axis
671 +  angle = dt * ji[2] / I[2][2];
672 +  this->rotate( 0, 1, angle, ji, A);
673 +  
674 +  // rotate about the y-axis
675 +  angle = dt2 * ji[1] / I[1][1];
676 +  this->rotate( 2, 0, angle, ji, A );
677 +  
678 +  // rotate about the x-axis
679 +  angle = dt2 * ji[0] / I[0][0];
680 +  this->rotate( 1, 2, angle, ji, A );
681 +  
682 +  dAtom->setA( A  );    
683   }
684  
685   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
# Line 749 | Line 772 | template<typename T> void Integrator<T>::thermalize(){
772   template<typename T> void Integrator<T>::thermalize(){
773    tStats->velocitize();
774   }
775 +
776 + template<typename T> double Integrator<T>::getConservedQuantity(void){
777 +  return tStats->getTotalE();
778 + }

Diff Legend

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