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 841 by mmeineke, Wed Oct 29 17:55:28 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 +  std::cerr << "Before initial Force calc\n";
178 +
179    calcForce(1, 1);
180 <  
180 >
181 >  if (nConstrained){
182 >    preMove();
183 >    constrainA();
184 >    calcForce(1, 1);
185 >    constrainB();
186 >    std::cerr << "premove done\n";
187 >  }
188 >
189 >
190 >
191    if (info->setTemp){
192      thermalize();
193    }
194  
182 –  calcPot = 0;
183 –  calcStress = 0;
184 –  currSample = sampleTime;
185 –  currThermal = thermalTime;
186 –  currStatus = statusTime;
187 –  
195    calcPot     = 0;
196    calcStress  = 0;
197    currSample  = sampleTime + info->getTime();
198    currThermal = thermalTime+ info->getTime();
199    currStatus  = statusTime + info->getTime();
200 +  currReset   = resetTime  + info->getTime();
201  
202    dumpOut->writeDump(info->getTime());
203    statOut->writeStat(info->getTime());
204  
197 –  readyCheck();
205  
206 +
207   #ifdef IS_MPI
208    strcpy(checkPointMsg, "The integrator is ready to go.");
209    MPIcheckPoint();
# Line 224 | Line 232 | template<typename T> void Integrator<T>::integrate(voi
232      }
233  
234      if (info->getTime() >= currStatus){
235 <      statOut->writeStat(info->getTime());
236 <      calcPot = 0;
235 >      statOut->writeStat(info->getTime());
236 >      calcPot = 0;
237        calcStress = 0;
238        currStatus += statusTime;
239 <    }
239 >    }
240  
241 +    if (info->resetIntegrator){
242 +      if (info->getTime() >= currReset){
243 +        this->resetIntegrator();
244 +        currReset += resetTime;
245 +      }
246 +    }
247 +
248 +    std::cerr << "done with time = " << info->getTime() << "\n";
249 +
250   #ifdef IS_MPI
251      strcpy(checkPointMsg, "successfully took a time step.");
252      MPIcheckPoint();
253   #endif // is_mpi
254    }
255  
239 –  dumpOut->writeFinal(info->getTime());
256  
257 +  // write the last frame
258 +  dumpOut->writeDump(info->getTime());
259 +
260    delete dumpOut;
261    delete statOut;
262   }
# Line 249 | Line 268 | template<typename T> void Integrator<T>::integrateStep
268  
269    moveA();
270  
252 –  if (nConstrained){
253 –    constrainA();
254 –  }
271  
272  
273 +
274   #ifdef IS_MPI
275    strcpy(checkPointMsg, "Succesful moveA\n");
276    MPIcheckPoint();
# Line 274 | Line 291 | template<typename T> void Integrator<T>::integrateStep
291  
292    moveB();
293  
277 –  if (nConstrained){
278 –    constrainB();
279 –  }
294  
295 +
296   #ifdef IS_MPI
297    strcpy(checkPointMsg, "Succesful moveB\n");
298    MPIcheckPoint();
# Line 289 | Line 304 | template<typename T> void Integrator<T>::moveA(void){
304    int i, j;
305    DirectionalAtom* dAtom;
306    double Tb[3], ji[3];
292 –  double A[3][3], I[3][3];
293 –  double angle;
307    double vel[3], pos[3], frc[3];
308    double mass;
309  
# Line 326 | Line 339 | template<typename T> void Integrator<T>::moveA(void){
339        for (j = 0; j < 3; j++)
340          ji[j] += (dt2 * Tb[j]) * eConvert;
341  
342 <      // use the angular velocities to propagate the rotation matrix a
330 <      // full time step
342 >      this->rotationPropagation( dAtom, ji );
343  
344 <      dAtom->getA(A);
345 <      dAtom->getI(I);
346 <
335 <      // rotate about the x-axis      
336 <      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);
344 >      dAtom->setJ(ji);
345 >    }
346 >  }
347  
348 <      // rotate about the z-axis
349 <      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 <    }
348 >  if (nConstrained){
349 >    constrainA();
350    }
351   }
352  
# Line 382 | Line 373 | template<typename T> void Integrator<T>::moveB(void){
373      if (atoms[i]->isDirectional()){
374        dAtom = (DirectionalAtom *) atoms[i];
375  
376 <      // get and convert the torque to body frame      
376 >      // get and convert the torque to body frame
377  
378        dAtom->getTrq(Tb);
379        dAtom->lab2Body(Tb);
# Line 398 | Line 389 | template<typename T> void Integrator<T>::moveB(void){
389        dAtom->setJ(ji);
390      }
391    }
392 +
393 +  if (nConstrained){
394 +    constrainB();
395 +  }
396   }
397  
398   template<typename T> void Integrator<T>::preMove(void){
# Line 416 | Line 411 | template<typename T> void Integrator<T>::constrainA(){
411   }
412  
413   template<typename T> void Integrator<T>::constrainA(){
414 <  int i, j, k;
414 >  int i, j;
415    int done;
416    double posA[3], posB[3];
417    double velA[3], velB[3];
# Line 556 | Line 551 | template<typename T> void Integrator<T>::constrainA(){
551      painCave.isFatal = 1;
552      simError();
553    }
554 +
555   }
556  
557   template<typename T> void Integrator<T>::constrainB(void){
558 <  int i, j, k;
558 >  int i, j;
559    int done;
560    double posA[3], posB[3];
561    double velA[3], velB[3];
# Line 568 | Line 564 | template<typename T> void Integrator<T>::constrainB(vo
564    int a, b, ax, ay, az, bx, by, bz;
565    double rma, rmb;
566    double dx, dy, dz;
567 <  double rabsq, pabsq, rvab;
572 <  double diffsq;
567 >  double rvab;
568    double gab;
569    int iteration;
570  
# Line 659 | Line 654 | template<typename T> void Integrator<T>::rotate(int ax
654    }
655   }
656  
657 + template<typename T> void Integrator<T>::rotationPropagation
658 + ( DirectionalAtom* dAtom, double ji[3] ){
659 +
660 +  double angle;
661 +  double A[3][3], I[3][3];
662 +
663 +  // use the angular velocities to propagate the rotation matrix a
664 +  // full time step
665 +
666 +  dAtom->getA(A);
667 +  dAtom->getI(I);
668 +
669 +  // rotate about the x-axis
670 +  angle = dt2 * ji[0] / I[0][0];
671 +  this->rotate( 1, 2, 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 z-axis
678 +  angle = dt * ji[2] / I[2][2];
679 +  this->rotate( 0, 1, angle, ji, A);
680 +
681 +  // rotate about the y-axis
682 +  angle = dt2 * ji[1] / I[1][1];
683 +  this->rotate( 2, 0, angle, ji, A );
684 +
685 +  // rotate about the x-axis
686 +  angle = dt2 * ji[0] / I[0][0];
687 +  this->rotate( 1, 2, angle, ji, A );
688 +
689 +  dAtom->setA( A  );
690 + }
691 +
692   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
693                                                  double angle, double ji[3],
694                                                  double A[3][3]){
# Line 724 | Line 754 | template<typename T> void Integrator<T>::rotate(int ax
754      }
755    }
756  
757 <  // rotate the Rotation matrix acording to:
757 >  // rotate the Rotation matrix acording to:
758    //            A[][] = A[][] * transpose(rot[][])
759  
760  
# Line 749 | Line 779 | template<typename T> void Integrator<T>::thermalize(){
779   template<typename T> void Integrator<T>::thermalize(){
780    tStats->velocitize();
781   }
782 +
783 + template<typename T> double Integrator<T>::getConservedQuantity(void){
784 +  return tStats->getTotalE();
785 + }
786 + template<typename T> string Integrator<T>::getAdditionalParameters(void){
787 +  //By default, return a null string
788 +  //The reason we use string instead of char* is that if we use char*, we will
789 +  //return a pointer point to local variable which might cause problem
790 +  return string();
791 + }

Diff Legend

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