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 787 by mmeineke, Thu Sep 25 19:27:15 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 147 | Line 146 | template<typename T> void Integrator<T>::integrate(voi
146  
147  
148   template<typename T> void Integrator<T>::integrate(void){
150 –  int i, j;                         // loop counters
149  
150    double runTime = info->run_time;
151    double sampleTime = info->sampleTime;
152    double statusTime = info->statusTime;
153    double thermalTime = info->thermalTime;
154 +  double resetTime = info->resetTime;
155  
156 +
157    double currSample;
158    double currThermal;
159    double currStatus;
160 <
160 >  double currReset;
161 >  
162    int calcPot, calcStress;
162 –  int isError;
163  
164    tStats = new Thermo(info);
165    statOut = new StatWriter(info);
166    dumpOut = new DumpWriter(info);
167  
168    atoms = info->atoms;
169 –  DirectionalAtom* dAtom;
169  
170    dt = info->dt;
171    dt2 = 0.5 * dt;
172  
173 +  readyCheck();
174 +
175    // initialize the forces before the first step
176  
177    calcForce(1, 1);
177 –  // myFF->doForces(1,1);
178  
179 +  if (nConstrained){
180 +    preMove();
181 +    constrainA();
182 +    calcForce(1, 1);    
183 +    constrainB();
184 +  }
185 +  
186    if (info->setTemp){
187      thermalize();
188    }
189  
183 –  calcPot = 0;
184 –  calcStress = 0;
185 –  currSample = sampleTime;
186 –  currThermal = thermalTime;
187 –  currStatus = statusTime;
188 –  
190    calcPot     = 0;
191    calcStress  = 0;
192    currSample  = sampleTime + info->getTime();
193    currThermal = thermalTime+ info->getTime();
194    currStatus  = statusTime + info->getTime();
195 +  currReset   = resetTime  + info->getTime();
196  
197    dumpOut->writeDump(info->getTime());
198    statOut->writeStat(info->getTime());
199  
198 –  readyCheck();
200  
201 +
202   #ifdef IS_MPI
203    strcpy(checkPointMsg, "The integrator is ready to go.");
204    MPIcheckPoint();
# Line 230 | Line 232 | template<typename T> void Integrator<T>::integrate(voi
232        calcStress = 0;
233        currStatus += statusTime;
234      }
235 +
236 +    if (info->resetIntegrator){
237 +      if (info->getTime() >= currReset){
238 +        this->resetIntegrator();
239 +        currReset += resetTime;
240 +      }
241 +    }
242  
243   #ifdef IS_MPI
244      strcpy(checkPointMsg, "successfully took a time step.");
# Line 250 | Line 259 | template<typename T> void Integrator<T>::integrateStep
259  
260    moveA();
261  
253 –  if (nConstrained){
254 –    constrainA();
255 –  }
262  
263  
264 +
265   #ifdef IS_MPI
266    strcpy(checkPointMsg, "Succesful moveA\n");
267    MPIcheckPoint();
# Line 275 | Line 282 | template<typename T> void Integrator<T>::integrateStep
282  
283    moveB();
284  
278 –  if (nConstrained){
279 –    constrainB();
280 –  }
285  
286 +
287   #ifdef IS_MPI
288    strcpy(checkPointMsg, "Succesful moveB\n");
289    MPIcheckPoint();
# Line 290 | Line 295 | template<typename T> void Integrator<T>::moveA(void){
295    int i, j;
296    DirectionalAtom* dAtom;
297    double Tb[3], ji[3];
293 –  double A[3][3], I[3][3];
294 –  double angle;
298    double vel[3], pos[3], frc[3];
299    double mass;
300  
# Line 327 | Line 330 | template<typename T> void Integrator<T>::moveA(void){
330        for (j = 0; j < 3; j++)
331          ji[j] += (dt2 * Tb[j]) * eConvert;
332  
333 <      // use the angular velocities to propagate the rotation matrix a
331 <      // full time step
333 >      this->rotationPropagation( dAtom, ji );
334  
335 <      dAtom->getA(A);
336 <      dAtom->getI(I);
337 <
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);
335 >      dAtom->setJ(ji);
336 >    }
337 >  }
338  
339 <      // rotate about the y-axis
340 <      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);
359 <    }
339 >  if (nConstrained){
340 >    constrainA();
341    }
342   }
343  
# Line 398 | Line 379 | template<typename T> void Integrator<T>::moveB(void){
379  
380        dAtom->setJ(ji);
381      }
382 +  }
383 +
384 +  if (nConstrained){
385 +    constrainB();
386    }
387   }
388  
# Line 417 | Line 402 | template<typename T> void Integrator<T>::constrainA(){
402   }
403  
404   template<typename T> void Integrator<T>::constrainA(){
405 <  int i, j, k;
405 >  int i, j;
406    int done;
407    double posA[3], posB[3];
408    double velA[3], velB[3];
# Line 557 | Line 542 | template<typename T> void Integrator<T>::constrainA(){
542      painCave.isFatal = 1;
543      simError();
544    }
545 +
546   }
547  
548   template<typename T> void Integrator<T>::constrainB(void){
549 <  int i, j, k;
549 >  int i, j;
550    int done;
551    double posA[3], posB[3];
552    double velA[3], velB[3];
# Line 569 | Line 555 | template<typename T> void Integrator<T>::constrainB(vo
555    int a, b, ax, ay, az, bx, by, bz;
556    double rma, rmb;
557    double dx, dy, dz;
558 <  double rabsq, pabsq, rvab;
573 <  double diffsq;
558 >  double rvab;
559    double gab;
560    int iteration;
561  
# Line 660 | Line 645 | template<typename T> void Integrator<T>::rotate(int ax
645    }
646   }
647  
648 + template<typename T> void Integrator<T>::rotationPropagation
649 + ( DirectionalAtom* dAtom, double ji[3] ){
650 +
651 +  double angle;
652 +  double A[3][3], I[3][3];
653 +
654 +  // use the angular velocities to propagate the rotation matrix a
655 +  // full time step
656 +
657 +  dAtom->getA(A);
658 +  dAtom->getI(I);
659 +  
660 +  // rotate about the x-axis      
661 +  angle = dt2 * ji[0] / I[0][0];
662 +  this->rotate( 1, 2, angle, ji, A );
663 +  
664 +  // rotate about the y-axis
665 +  angle = dt2 * ji[1] / I[1][1];
666 +  this->rotate( 2, 0, angle, ji, A );
667 +  
668 +  // rotate about the z-axis
669 +  angle = dt * ji[2] / I[2][2];
670 +  this->rotate( 0, 1, angle, ji, A);
671 +  
672 +  // rotate about the y-axis
673 +  angle = dt2 * ji[1] / I[1][1];
674 +  this->rotate( 2, 0, angle, ji, A );
675 +  
676 +  // rotate about the x-axis
677 +  angle = dt2 * ji[0] / I[0][0];
678 +  this->rotate( 1, 2, angle, ji, A );
679 +  
680 +  dAtom->setA( A  );    
681 + }
682 +
683   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
684                                                  double angle, double ji[3],
685                                                  double A[3][3]){
# Line 750 | Line 770 | template<typename T> void Integrator<T>::thermalize(){
770   template<typename T> void Integrator<T>::thermalize(){
771    tStats->velocitize();
772   }
773 +
774 + template<typename T> double Integrator<T>::getConservedQuantity(void){
775 +  return tStats->getTotalE();
776 + }

Diff Legend

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