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 829 by gezelter, Tue Oct 28 16:03:37 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;
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);
178 +
179 +  if (nConstrained){
180 +    preMove();
181 +    constrainA();
182 +    calcForce(1, 1);    
183 +    constrainB();
184 +  }
185    
186    if (info->setTemp){
187      thermalize();
188    }
189  
182 –  calcPot = 0;
183 –  calcStress = 0;
184 –  currSample = sampleTime;
185 –  currThermal = thermalTime;
186 –  currStatus = statusTime;
187 –  
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  
197 –  readyCheck();
200  
201 +
202   #ifdef IS_MPI
203    strcpy(checkPointMsg, "The integrator is ready to go.");
204    MPIcheckPoint();
# Line 230 | Line 233 | template<typename T> void Integrator<T>::integrate(voi
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.");
245      MPIcheckPoint();
246   #endif // is_mpi
247    }
248  
239 –  dumpOut->writeFinal(info->getTime());
249  
250 +  // write the last frame
251 +  dumpOut->writeDump(info->getTime());
252 +
253    delete dumpOut;
254    delete statOut;
255   }
# Line 249 | Line 261 | template<typename T> void Integrator<T>::integrateStep
261  
262    moveA();
263  
252 –  if (nConstrained){
253 –    constrainA();
254 –  }
264  
265  
266 +
267   #ifdef IS_MPI
268    strcpy(checkPointMsg, "Succesful moveA\n");
269    MPIcheckPoint();
# Line 274 | Line 284 | template<typename T> void Integrator<T>::integrateStep
284  
285    moveB();
286  
277 –  if (nConstrained){
278 –    constrainB();
279 –  }
287  
288 +
289   #ifdef IS_MPI
290    strcpy(checkPointMsg, "Succesful moveB\n");
291    MPIcheckPoint();
# Line 289 | Line 297 | template<typename T> void Integrator<T>::moveA(void){
297    int i, j;
298    DirectionalAtom* dAtom;
299    double Tb[3], ji[3];
292 –  double A[3][3], I[3][3];
293 –  double angle;
300    double vel[3], pos[3], frc[3];
301    double mass;
302  
# Line 326 | 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
330 <      // full time step
335 >      this->rotationPropagation( dAtom, ji );
336  
337 <      dAtom->getA(A);
338 <      dAtom->getI(I);
337 >      dAtom->setJ(ji);
338 >    }
339 >  }
340  
341 <      // rotate about the x-axis      
342 <      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);
342 <
343 <      // rotate about the z-axis
344 <      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 <    }
341 >  if (nConstrained){
342 >    constrainA();
343    }
344   }
345  
# Line 398 | 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 416 | Line 404 | template<typename T> void Integrator<T>::constrainA(){
404   }
405  
406   template<typename T> void Integrator<T>::constrainA(){
407 <  int i, j, k;
407 >  int i, j;
408    int done;
409    double posA[3], posB[3];
410    double velA[3], velB[3];
# Line 556 | 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){
551 <  int i, j, k;
551 >  int i, j;
552    int done;
553    double posA[3], posB[3];
554    double velA[3], velB[3];
# Line 568 | Line 557 | template<typename T> void Integrator<T>::constrainB(vo
557    int a, b, ax, ay, az, bx, by, bz;
558    double rma, rmb;
559    double dx, dy, dz;
560 <  double rabsq, pabsq, rvab;
572 <  double diffsq;
560 >  double rvab;
561    double gab;
562    int iteration;
563  
# Line 659 | Line 647 | template<typename T> void Integrator<T>::rotate(int ax
647    }
648   }
649  
650 + template<typename T> void Integrator<T>::rotationPropagation
651 + ( DirectionalAtom* dAtom, double ji[3] ){
652 +
653 +  double angle;
654 +  double A[3][3], I[3][3];
655 +
656 +  // use the angular velocities to propagate the rotation matrix a
657 +  // full time step
658 +
659 +  dAtom->getA(A);
660 +  dAtom->getI(I);
661 +  
662 +  // rotate about the x-axis      
663 +  angle = dt2 * ji[0] / I[0][0];
664 +  this->rotate( 1, 2, angle, ji, A );
665 +  
666 +  // rotate about the y-axis
667 +  angle = dt2 * ji[1] / I[1][1];
668 +  this->rotate( 2, 0, angle, ji, A );
669 +  
670 +  // rotate about the z-axis
671 +  angle = dt * ji[2] / I[2][2];
672 +  this->rotate( 0, 1, angle, ji, A);
673 +  
674 +  // rotate about the y-axis
675 +  angle = dt2 * ji[1] / I[1][1];
676 +  this->rotate( 2, 0, angle, ji, A );
677 +  
678 +  // rotate about the x-axis
679 +  angle = dt2 * ji[0] / I[0][0];
680 +  this->rotate( 1, 2, angle, ji, A );
681 +  
682 +  dAtom->setA( A  );    
683 + }
684 +
685   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
686                                                  double angle, double ji[3],
687                                                  double A[3][3]){
# Line 749 | Line 772 | template<typename T> void Integrator<T>::thermalize(){
772   template<typename T> void Integrator<T>::thermalize(){
773    tStats->velocitize();
774   }
775 +
776 + template<typename T> double Integrator<T>::getConservedQuantity(void){
777 +  return tStats->getTotalE();
778 + }

Diff Legend

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