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 892 by chuckv, Mon Dec 22 21:27:04 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"
7   #include <unistd.h>
8   #endif //is_mpi
9  
10 + #ifdef PROFILE
11 + #include "mdProfile.hpp"
12 + #endif // profile
13 +
14   #include "Integrator.hpp"
15   #include "simError.h"
16  
# Line 25 | Line 29 | template<typename T> Integrator<T>::Integrator(SimInfo
29    if (info->the_integrator != NULL){
30      delete info->the_integrator;
31    }
28 –  info->the_integrator = this;
32  
33    nAtoms = info->n_atoms;
34  
# Line 147 | Line 150 | template<typename T> void Integrator<T>::integrate(voi
150  
151  
152   template<typename T> void Integrator<T>::integrate(void){
150 –  int i, j;                         // loop counters
153  
154    double runTime = info->run_time;
155    double sampleTime = info->sampleTime;
156    double statusTime = info->statusTime;
157    double thermalTime = info->thermalTime;
158 +  double resetTime = info->resetTime;
159  
160 +
161    double currSample;
162    double currThermal;
163    double currStatus;
164 +  double currReset;
165  
166    int calcPot, calcStress;
162 –  int isError;
167  
168    tStats = new Thermo(info);
169    statOut = new StatWriter(info);
170    dumpOut = new DumpWriter(info);
171  
172    atoms = info->atoms;
169 –  DirectionalAtom* dAtom;
173  
174    dt = info->dt;
175    dt2 = 0.5 * dt;
176  
177 +  readyCheck();
178 +
179    // initialize the forces before the first step
180  
181    calcForce(1, 1);
177 –  // myFF->doForces(1,1);
182  
183 +  if (nConstrained){
184 +    preMove();
185 +    constrainA();
186 +    calcForce(1, 1);
187 +    constrainB();
188 +  }
189 +  
190    if (info->setTemp){
191      thermalize();
192    }
193  
183 –  calcPot = 0;
184 –  calcStress = 0;
185 –  currSample = sampleTime;
186 –  currThermal = thermalTime;
187 –  currStatus = statusTime;
188 –  
194    calcPot     = 0;
195    calcStress  = 0;
196    currSample  = sampleTime + info->getTime();
197    currThermal = thermalTime+ info->getTime();
198    currStatus  = statusTime + info->getTime();
199 +  currReset   = resetTime  + info->getTime();
200  
201    dumpOut->writeDump(info->getTime());
202    statOut->writeStat(info->getTime());
203  
198 –  readyCheck();
204  
205   #ifdef IS_MPI
206    strcpy(checkPointMsg, "The integrator is ready to go.");
# Line 208 | Line 213 | template<typename T> void Integrator<T>::integrate(voi
213        calcStress = 1;
214      }
215  
216 + #ifdef PROFILE
217 +    startProfile( pro1 );
218 + #endif
219 +    
220      integrateStep(calcPot, calcStress);
221  
222 + #ifdef PROFILE
223 +    endProfile( pro1 );
224 +
225 +    startProfile( pro2 );
226 + #endif // profile
227 +
228      info->incrTime(dt);
229  
230      if (info->setTemp){
# Line 225 | Line 240 | template<typename T> void Integrator<T>::integrate(voi
240      }
241  
242      if (info->getTime() >= currStatus){
243 <      statOut->writeStat(info->getTime());
244 <      calcPot = 0;
243 >      statOut->writeStat(info->getTime());
244 >      calcPot = 0;
245        calcStress = 0;
246        currStatus += statusTime;
247 <    }
247 >    }
248  
249 +    if (info->resetIntegrator){
250 +      if (info->getTime() >= currReset){
251 +        this->resetIntegrator();
252 +        currReset += resetTime;
253 +      }
254 +    }
255 +    
256 + #ifdef PROFILE
257 +    endProfile( pro2 );
258 + #endif //profile
259 +
260   #ifdef IS_MPI
261      strcpy(checkPointMsg, "successfully took a time step.");
262      MPIcheckPoint();
263   #endif // is_mpi
264    }
265  
240 –  dumpOut->writeFinal(info->getTime());
266  
267 +  // write the last frame
268 +  dumpOut->writeDump(info->getTime());
269 +
270    delete dumpOut;
271    delete statOut;
272   }
# Line 246 | Line 274 | template<typename T> void Integrator<T>::integrateStep
274   template<typename T> void Integrator<T>::integrateStep(int calcPot,
275                                                         int calcStress){
276    // Position full step, and velocity half step
277 +
278 + #ifdef PROFILE
279 +  startProfile(pro3);
280 + #endif //profile
281 +
282    preMove();
283  
284 + #ifdef PROFILE
285 +  endProfile(pro3);
286 +
287 +  startProfile(pro4);
288 + #endif // profile
289 +
290    moveA();
291  
292 <  if (nConstrained){
293 <    constrainA();
294 <  }
292 > #ifdef PROFILE
293 >  endProfile(pro4);
294 >  
295 >  startProfile(pro5);
296 > #endif//profile
297  
298  
299   #ifdef IS_MPI
# Line 270 | Line 311 | template<typename T> void Integrator<T>::integrateStep
311    MPIcheckPoint();
312   #endif // is_mpi
313  
314 + #ifdef PROFILE
315 +  endProfile( pro5 );
316  
317 +  startProfile( pro6 );
318 + #endif //profile
319 +
320    // finish the velocity  half step
321  
322    moveB();
323  
324 <  if (nConstrained){
325 <    constrainB();
326 <  }
324 > #ifdef PROFILE
325 >  endProfile(pro6);
326 > #endif // profile
327  
328   #ifdef IS_MPI
329    strcpy(checkPointMsg, "Succesful moveB\n");
# Line 290 | Line 336 | template<typename T> void Integrator<T>::moveA(void){
336    int i, j;
337    DirectionalAtom* dAtom;
338    double Tb[3], ji[3];
293 –  double A[3][3], I[3][3];
294 –  double angle;
339    double vel[3], pos[3], frc[3];
340    double mass;
341  
# Line 327 | Line 371 | template<typename T> void Integrator<T>::moveA(void){
371        for (j = 0; j < 3; j++)
372          ji[j] += (dt2 * Tb[j]) * eConvert;
373  
374 <      // use the angular velocities to propagate the rotation matrix a
331 <      // full time step
332 <
333 <      dAtom->getA(A);
334 <      dAtom->getI(I);
335 <
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);
347 <
348 <      // rotate about the y-axis
349 <      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 <
374 >      this->rotationPropagation( dAtom, ji );
375  
376        dAtom->setJ(ji);
358 –      dAtom->setA(A);
377      }
378    }
379 +
380 +  if (nConstrained){
381 +    constrainA();
382 +  }
383   }
384  
385  
# Line 383 | Line 405 | template<typename T> void Integrator<T>::moveB(void){
405      if (atoms[i]->isDirectional()){
406        dAtom = (DirectionalAtom *) atoms[i];
407  
408 <      // get and convert the torque to body frame      
408 >      // get and convert the torque to body frame
409  
410        dAtom->getTrq(Tb);
411        dAtom->lab2Body(Tb);
# Line 399 | Line 421 | template<typename T> void Integrator<T>::moveB(void){
421        dAtom->setJ(ji);
422      }
423    }
424 +
425 +  if (nConstrained){
426 +    constrainB();
427 +  }
428   }
429  
430   template<typename T> void Integrator<T>::preMove(void){
# Line 417 | Line 443 | template<typename T> void Integrator<T>::constrainA(){
443   }
444  
445   template<typename T> void Integrator<T>::constrainA(){
446 <  int i, j, k;
446 >  int i, j;
447    int done;
448    double posA[3], posB[3];
449    double velA[3], velB[3];
# Line 557 | Line 583 | template<typename T> void Integrator<T>::constrainA(){
583      painCave.isFatal = 1;
584      simError();
585    }
586 +
587   }
588  
589   template<typename T> void Integrator<T>::constrainB(void){
590 <  int i, j, k;
590 >  int i, j;
591    int done;
592    double posA[3], posB[3];
593    double velA[3], velB[3];
# Line 569 | Line 596 | template<typename T> void Integrator<T>::constrainB(vo
596    int a, b, ax, ay, az, bx, by, bz;
597    double rma, rmb;
598    double dx, dy, dz;
599 <  double rabsq, pabsq, rvab;
573 <  double diffsq;
599 >  double rvab;
600    double gab;
601    int iteration;
602  
# Line 660 | Line 686 | template<typename T> void Integrator<T>::rotate(int ax
686    }
687   }
688  
689 + template<typename T> void Integrator<T>::rotationPropagation
690 + ( DirectionalAtom* dAtom, double ji[3] ){
691 +
692 +  double angle;
693 +  double A[3][3], I[3][3];
694 +
695 +  // use the angular velocities to propagate the rotation matrix a
696 +  // full time step
697 +
698 +  dAtom->getA(A);
699 +  dAtom->getI(I);
700 +
701 +  // rotate about the x-axis
702 +  angle = dt2 * ji[0] / I[0][0];
703 +  this->rotate( 1, 2, angle, ji, A );
704 +
705 +  // rotate about the y-axis
706 +  angle = dt2 * ji[1] / I[1][1];
707 +  this->rotate( 2, 0, angle, ji, A );
708 +
709 +  // rotate about the z-axis
710 +  angle = dt * ji[2] / I[2][2];
711 +  this->rotate( 0, 1, angle, ji, A);
712 +
713 +  // rotate about the y-axis
714 +  angle = dt2 * ji[1] / I[1][1];
715 +  this->rotate( 2, 0, angle, ji, A );
716 +
717 +  // rotate about the x-axis
718 +  angle = dt2 * ji[0] / I[0][0];
719 +  this->rotate( 1, 2, angle, ji, A );
720 +
721 +  dAtom->setA( A  );
722 + }
723 +
724   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
725                                                  double angle, double ji[3],
726                                                  double A[3][3]){
# Line 725 | Line 786 | template<typename T> void Integrator<T>::rotate(int ax
786      }
787    }
788  
789 <  // rotate the Rotation matrix acording to:
789 >  // rotate the Rotation matrix acording to:
790    //            A[][] = A[][] * transpose(rot[][])
791  
792  
# Line 750 | Line 811 | template<typename T> void Integrator<T>::thermalize(){
811   template<typename T> void Integrator<T>::thermalize(){
812    tStats->velocitize();
813   }
814 +
815 + template<typename T> double Integrator<T>::getConservedQuantity(void){
816 +  return tStats->getTotalE();
817 + }
818 + template<typename T> string Integrator<T>::getAdditionalParameters(void){
819 +  //By default, return a null string
820 +  //The reason we use string instead of char* is that if we use char*, we will
821 +  //return a pointer point to local variable which might cause problem
822 +  return string();
823 + }

Diff Legend

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