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 843 by mmeineke, Wed Oct 29 20:41:39 2003 UTC

# Line 1 | Line 1
1   #include <iostream>
2 < #include <cstdlib>
3 < #include <cmath>
2 > #include <stdlib.h>
3 > #include <math.h>
4  
5   #ifdef IS_MPI
6   #include "mpiSimulation.hpp"
# 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;
28  
29    nAtoms = info->n_atoms;
30  
# 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 +  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 < >>>>>>> 1.18
195 >  currReset   = resetTime  + info->getTime();
196  
197    dumpOut->writeDump(info->getTime());
198    statOut->writeStat(info->getTime());
199  
199 –  readyCheck();
200  
201   #ifdef IS_MPI
202    strcpy(checkPointMsg, "The integrator is ready to go.");
# Line 226 | Line 226 | template<typename T> void Integrator<T>::integrate(voi
226      }
227  
228      if (info->getTime() >= currStatus){
229 <      statOut->writeStat(info->getTime());
230 <      calcPot = 0;
229 >      statOut->writeStat(info->getTime());
230 >      calcPot = 0;
231        calcStress = 0;
232        currStatus += statusTime;
233 <    }
233 >    }
234 >
235 >    if (info->resetIntegrator){
236 >      if (info->getTime() >= currReset){
237 >        this->resetIntegrator();
238 >        currReset += resetTime;
239 >      }
240 >    }
241  
242   #ifdef IS_MPI
243      strcpy(checkPointMsg, "successfully took a time step.");
# Line 238 | Line 245 | template<typename T> void Integrator<T>::integrate(voi
245   #endif // is_mpi
246    }
247  
241 –  dumpOut->writeFinal(info->getTime());
248  
249 +  // write the last frame
250 +  dumpOut->writeDump(info->getTime());
251 +
252    delete dumpOut;
253    delete statOut;
254   }
# Line 251 | Line 260 | template<typename T> void Integrator<T>::integrateStep
260  
261    moveA();
262  
254 –  if (nConstrained){
255 –    constrainA();
256 –  }
263  
264  
265 +
266   #ifdef IS_MPI
267    strcpy(checkPointMsg, "Succesful moveA\n");
268    MPIcheckPoint();
# Line 276 | Line 283 | template<typename T> void Integrator<T>::integrateStep
283  
284    moveB();
285  
279 –  if (nConstrained){
280 –    constrainB();
281 –  }
286  
287 +
288   #ifdef IS_MPI
289    strcpy(checkPointMsg, "Succesful moveB\n");
290    MPIcheckPoint();
# Line 291 | Line 296 | template<typename T> void Integrator<T>::moveA(void){
296    int i, j;
297    DirectionalAtom* dAtom;
298    double Tb[3], ji[3];
294 –  double A[3][3], I[3][3];
295 –  double angle;
299    double vel[3], pos[3], frc[3];
300    double mass;
301  
# Line 328 | 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
332 <      // full time step
334 >      this->rotationPropagation( dAtom, ji );
335  
336 <      dAtom->getA(A);
337 <      dAtom->getI(I);
338 <
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);
336 >      dAtom->setJ(ji);
337 >    }
338 >  }
339  
340 <      // rotate about the y-axis
341 <      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 <    }
340 >  if (nConstrained){
341 >    constrainA();
342    }
343   }
344  
# Line 384 | Line 365 | template<typename T> void Integrator<T>::moveB(void){
365      if (atoms[i]->isDirectional()){
366        dAtom = (DirectionalAtom *) atoms[i];
367  
368 <      // get and convert the torque to body frame      
368 >      // get and convert the torque to body frame
369  
370        dAtom->getTrq(Tb);
371        dAtom->lab2Body(Tb);
# Line 400 | 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 418 | Line 403 | template<typename T> void Integrator<T>::constrainA(){
403   }
404  
405   template<typename T> void Integrator<T>::constrainA(){
406 <  int i, j, k;
406 >  int i, j;
407    int done;
408    double posA[3], posB[3];
409    double velA[3], velB[3];
# Line 558 | 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){
550 <  int i, j, k;
550 >  int i, j;
551    int done;
552    double posA[3], posB[3];
553    double velA[3], velB[3];
# Line 570 | Line 556 | template<typename T> void Integrator<T>::constrainB(vo
556    int a, b, ax, ay, az, bx, by, bz;
557    double rma, rmb;
558    double dx, dy, dz;
559 <  double rabsq, pabsq, rvab;
574 <  double diffsq;
559 >  double rvab;
560    double gab;
561    int iteration;
562  
# Line 661 | Line 646 | template<typename T> void Integrator<T>::rotate(int ax
646    }
647   }
648  
649 + template<typename T> void Integrator<T>::rotationPropagation
650 + ( DirectionalAtom* dAtom, double ji[3] ){
651 +
652 +  double angle;
653 +  double A[3][3], I[3][3];
654 +
655 +  // use the angular velocities to propagate the rotation matrix a
656 +  // full time step
657 +
658 +  dAtom->getA(A);
659 +  dAtom->getI(I);
660 +
661 +  // rotate about the x-axis
662 +  angle = dt2 * ji[0] / I[0][0];
663 +  this->rotate( 1, 2, angle, ji, A );
664 +
665 +  // rotate about the y-axis
666 +  angle = dt2 * ji[1] / I[1][1];
667 +  this->rotate( 2, 0, angle, ji, A );
668 +
669 +  // rotate about the z-axis
670 +  angle = dt * ji[2] / I[2][2];
671 +  this->rotate( 0, 1, angle, ji, A);
672 +
673 +  // rotate about the y-axis
674 +  angle = dt2 * ji[1] / I[1][1];
675 +  this->rotate( 2, 0, angle, ji, A );
676 +
677 +  // rotate about the x-axis
678 +  angle = dt2 * ji[0] / I[0][0];
679 +  this->rotate( 1, 2, angle, ji, A );
680 +
681 +  dAtom->setA( A  );
682 + }
683 +
684   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
685                                                  double angle, double ji[3],
686                                                  double A[3][3]){
# Line 726 | Line 746 | template<typename T> void Integrator<T>::rotate(int ax
746      }
747    }
748  
749 <  // rotate the Rotation matrix acording to:
749 >  // rotate the Rotation matrix acording to:
750    //            A[][] = A[][] * transpose(rot[][])
751  
752  
# Line 751 | Line 771 | template<typename T> void Integrator<T>::thermalize(){
771   template<typename T> void Integrator<T>::thermalize(){
772    tStats->velocitize();
773   }
774 +
775 + template<typename T> double Integrator<T>::getConservedQuantity(void){
776 +  return tStats->getTotalE();
777 + }
778 + template<typename T> string Integrator<T>::getAdditionalParameters(void){
779 +  //By default, return a null string
780 +  //The reason we use string instead of char* is that if we use char*, we will
781 +  //return a pointer point to local variable which might cause problem
782 +  return string();
783 + }

Diff Legend

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