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 784 by mmeineke, Wed Sep 24 19:34:39 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 171 | Line 173 | template<typename T> void Integrator<T>::integrate(voi
173    dt = info->dt;
174    dt2 = 0.5 * dt;
175  
176 +  readyCheck();
177 +
178    // initialize the forces before the first step
179  
180    calcForce(1, 1);
177 –  // myFF->doForces(1,1);
181  
182 +  if (nConstrained){
183 +    preMove();
184 +    constrainA();
185 +    calcForce(1, 1);    
186 +    constrainB();
187 +  }
188 +  
189    if (info->setTemp){
190      thermalize();
191    }
192  
183 –  calcPot = 0;
184 –  calcStress = 0;
185 –  currSample = sampleTime;
186 –  currThermal = thermalTime;
187 –  currStatus = statusTime;
188 –  
193    calcPot     = 0;
194    calcStress  = 0;
195    currSample  = sampleTime + info->getTime();
196    currThermal = thermalTime+ info->getTime();
197    currStatus  = statusTime + info->getTime();
198 +  currReset   = resetTime  + info->getTime();
199  
200    dumpOut->writeDump(info->getTime());
201    statOut->writeStat(info->getTime());
202  
198 –  readyCheck();
203  
204 +
205   #ifdef IS_MPI
206    strcpy(checkPointMsg, "The integrator is ready to go.");
207    MPIcheckPoint();
# Line 231 | Line 236 | template<typename T> void Integrator<T>::integrate(voi
236        currStatus += statusTime;
237      }
238  
239 +    if (info->resetIntegrator){
240 +      if (info->getTime() >= currReset){
241 +        this->resetIntegrator();
242 +        currReset += resetTime;
243 +      }
244 +    }
245 +
246   #ifdef IS_MPI
247      strcpy(checkPointMsg, "successfully took a time step.");
248      MPIcheckPoint();
# Line 250 | Line 262 | template<typename T> void Integrator<T>::integrateStep
262  
263    moveA();
264  
253 –  if (nConstrained){
254 –    constrainA();
255 –  }
265  
266  
267 +
268   #ifdef IS_MPI
269    strcpy(checkPointMsg, "Succesful moveA\n");
270    MPIcheckPoint();
# Line 275 | Line 285 | template<typename T> void Integrator<T>::integrateStep
285  
286    moveB();
287  
278 –  if (nConstrained){
279 –    constrainB();
280 –  }
288  
289 +
290   #ifdef IS_MPI
291    strcpy(checkPointMsg, "Succesful moveB\n");
292    MPIcheckPoint();
# Line 290 | Line 298 | template<typename T> void Integrator<T>::moveA(void){
298    int i, j;
299    DirectionalAtom* dAtom;
300    double Tb[3], ji[3];
293 –  double A[3][3], I[3][3];
294 –  double angle;
301    double vel[3], pos[3], frc[3];
302    double mass;
303  
# Line 327 | Line 333 | template<typename T> void Integrator<T>::moveA(void){
333        for (j = 0; j < 3; j++)
334          ji[j] += (dt2 * Tb[j]) * eConvert;
335  
336 <      // use the angular velocities to propagate the rotation matrix a
331 <      // full time step
336 >      this->rotationPropagation( dAtom, ji );
337  
338 <      dAtom->getA(A);
334 <      dAtom->getI(I);
335 <
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);
351 <
352 <      // rotate about the x-axis
353 <      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);
338 >      dAtom->setJ(ji);
339      }
340 +  }
341 +
342 +  if (nConstrained){
343 +    constrainA();
344    }
345   }
346  
# Line 399 | Line 383 | template<typename T> void Integrator<T>::moveB(void){
383        dAtom->setJ(ji);
384      }
385    }
386 +
387 +  if (nConstrained){
388 +    constrainB();
389 +  }
390   }
391  
392   template<typename T> void Integrator<T>::preMove(void){
# Line 557 | Line 545 | template<typename T> void Integrator<T>::constrainA(){
545      painCave.isFatal = 1;
546      simError();
547    }
548 +
549   }
550  
551   template<typename T> void Integrator<T>::constrainB(void){
# Line 658 | Line 647 | template<typename T> void Integrator<T>::constrainB(vo
647      painCave.isFatal = 1;
648      simError();
649    }
650 + }
651 +
652 + template<typename T> void Integrator<T>::rotationPropagation
653 + ( DirectionalAtom* dAtom, double ji[3] ){
654 +
655 +  double angle;
656 +  double A[3][3], I[3][3];
657 +
658 +  // use the angular velocities to propagate the rotation matrix a
659 +  // full time step
660 +
661 +  dAtom->getA(A);
662 +  dAtom->getI(I);
663 +  
664 +  // rotate about the x-axis      
665 +  angle = dt2 * ji[0] / I[0][0];
666 +  this->rotate( 1, 2, 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 z-axis
673 +  angle = dt * ji[2] / I[2][2];
674 +  this->rotate( 0, 1, angle, ji, A);
675 +  
676 +  // rotate about the y-axis
677 +  angle = dt2 * ji[1] / I[1][1];
678 +  this->rotate( 2, 0, angle, ji, A );
679 +  
680 +  // rotate about the x-axis
681 +  angle = dt2 * ji[0] / I[0][0];
682 +  this->rotate( 1, 2, angle, ji, A );
683 +  
684 +  dAtom->setA( A  );    
685   }
686  
687   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
# Line 750 | Line 774 | template<typename T> void Integrator<T>::thermalize(){
774   template<typename T> void Integrator<T>::thermalize(){
775    tStats->velocitize();
776   }
777 +
778 + template<typename T> double Integrator<T>::getConservedQuantity(void){
779 +  return tStats->getTotalE();
780 + }

Diff Legend

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