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 725 by tim, Tue Aug 26 20:29:26 2003 UTC vs.
Revision 781 by tim, Mon Sep 22 23:07:57 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);
177 –  // myFF->doForces(1,1);
180  
181 +  if (nConstrained){
182 +    preMove();
183 +    constrainA();
184 +    calcForce(1, 1);    
185 +    constrainB();
186 +  }
187 +  
188    if (info->setTemp){
189      thermalize();
190    }
191  
183 –  calcPot = 0;
184 –  calcStress = 0;
185 –  currSample = sampleTime;
186 –  currThermal = thermalTime;
187 –  currStatus = statusTime;
188 –  
192    calcPot     = 0;
193    calcStress  = 0;
194    currSample  = sampleTime + info->getTime();
195    currThermal = thermalTime+ info->getTime();
196    currStatus  = statusTime + info->getTime();
197 < >>>>>>> 1.18
197 >  currReset   = resetTime  + info->getTime();
198  
199    dumpOut->writeDump(info->getTime());
200    statOut->writeStat(info->getTime());
# Line 232 | Line 235 | template<typename T> void Integrator<T>::integrate(voi
235        currStatus += statusTime;
236      }
237  
238 +    if (info->resetIntegrator){
239 +      if (info->getTime() >= currReset){
240 +        this->resetIntegrator();
241 +        currReset += resetTime;
242 +      }
243 +    }
244 +
245   #ifdef IS_MPI
246      strcpy(checkPointMsg, "successfully took a time step.");
247      MPIcheckPoint();
# Line 251 | Line 261 | template<typename T> void Integrator<T>::integrateStep
261  
262    moveA();
263  
264 <  if (nConstrained){
255 <    constrainA();
256 <  }
264 >
265  
266  
267   #ifdef IS_MPI
# Line 276 | Line 284 | template<typename T> void Integrator<T>::integrateStep
284  
285    moveB();
286  
279 –  if (nConstrained){
280 –    constrainB();
281 –  }
287  
288 +
289   #ifdef IS_MPI
290    strcpy(checkPointMsg, "Succesful moveB\n");
291    MPIcheckPoint();
# Line 291 | Line 297 | template<typename T> void Integrator<T>::moveA(void){
297    int i, j;
298    DirectionalAtom* dAtom;
299    double Tb[3], ji[3];
294 –  double A[3][3], I[3][3];
295 –  double angle;
300    double vel[3], pos[3], frc[3];
301    double mass;
302  
# Line 328 | Line 332 | template<typename T> void Integrator<T>::moveA(void){
332        for (j = 0; j < 3; j++)
333          ji[j] += (dt2 * Tb[j]) * eConvert;
334  
335 <      // use the angular velocities to propagate the rotation matrix a
332 <      // full time step
335 >      this->rotationPropagation( dAtom, ji );
336  
337 <      dAtom->getA(A);
338 <      dAtom->getI(I);
339 <
337 <      // rotate about the x-axis      
338 <      angle = dt2 * ji[0] / I[0][0];
339 <      this->rotate(1, 2, angle, ji, A);
340 <
341 <      // rotate about the y-axis
342 <      angle = dt2 * ji[1] / I[1][1];
343 <      this->rotate(2, 0, angle, ji, A);
344 <
345 <      // rotate about the z-axis
346 <      angle = dt * ji[2] / I[2][2];
347 <      this->rotate(0, 1, angle, ji, A);
337 >      dAtom->setJ(ji);
338 >    }
339 >  }
340  
341 <      // rotate about the y-axis
342 <      angle = dt2 * ji[1] / I[1][1];
351 <      this->rotate(2, 0, angle, ji, A);
352 <
353 <      // rotate about the x-axis
354 <      angle = dt2 * ji[0] / I[0][0];
355 <      this->rotate(1, 2, angle, ji, A);
356 <
357 <
358 <      dAtom->setJ(ji);
359 <      dAtom->setA(A);
360 <    }
341 >  if (nConstrained){
342 >    constrainA();
343    }
344   }
345  
# Line 400 | Line 382 | template<typename T> void Integrator<T>::moveB(void){
382        dAtom->setJ(ji);
383      }
384    }
385 +
386 +  if (nConstrained){
387 +    constrainB();
388 +  }
389   }
390  
391   template<typename T> void Integrator<T>::preMove(void){
# Line 558 | Line 544 | template<typename T> void Integrator<T>::constrainA(){
544      painCave.isFatal = 1;
545      simError();
546    }
547 +
548   }
549  
550   template<typename T> void Integrator<T>::constrainB(void){
# Line 659 | Line 646 | template<typename T> void Integrator<T>::constrainB(vo
646      painCave.isFatal = 1;
647      simError();
648    }
649 + }
650 +
651 + template<typename T> void Integrator<T>::rotationPropagation
652 + ( DirectionalAtom* dAtom, double ji[3] ){
653 +
654 +  double angle;
655 +  double A[3][3], I[3][3];
656 +
657 +  // use the angular velocities to propagate the rotation matrix a
658 +  // full time step
659 +
660 +  dAtom->getA(A);
661 +  dAtom->getI(I);
662 +  
663 +  // rotate about the x-axis      
664 +  angle = dt2 * ji[0] / I[0][0];
665 +  this->rotate( 1, 2, angle, ji, A );
666 +  
667 +  // rotate about the y-axis
668 +  angle = dt2 * ji[1] / I[1][1];
669 +  this->rotate( 2, 0, angle, ji, A );
670 +  
671 +  // rotate about the z-axis
672 +  angle = dt * ji[2] / I[2][2];
673 +  this->rotate( 0, 1, angle, ji, A);
674 +  
675 +  // rotate about the y-axis
676 +  angle = dt2 * ji[1] / I[1][1];
677 +  this->rotate( 2, 0, angle, ji, A );
678 +  
679 +  // rotate about the x-axis
680 +  angle = dt2 * ji[0] / I[0][0];
681 +  this->rotate( 1, 2, angle, ji, A );
682 +  
683 +  dAtom->setA( A  );    
684   }
685  
686   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
# Line 751 | Line 773 | template<typename T> void Integrator<T>::thermalize(){
773   template<typename T> void Integrator<T>::thermalize(){
774    tStats->velocitize();
775   }
776 +
777 + template<typename T> double Integrator<T>::getConservedQuantity(void){
778 +  return tStats->getTotalE();
779 + }

Diff Legend

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